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1  Summary 

This  report  summarizes  the  work  completed  during  the  third  year  of  the  three  year  AFOSR 
Grant  no.  F  49620-03-1-0298  to  the  University  of  Kansas,  Lawrence,  KS  (K.S.Surana, 
PI)  and  Grant  no.  F  49620-03-1-0201  to  Texas  A  k,  M  University,  College  Station,  TX 
(J.N. Reddy,  PI).  The  development  of  hpk  mathematical  and  computational  framework 
for  all  boundary  value  problems  (BVPs)  and  initial  value  problems  (IVPs)  regardless  of 
their  origin  or  fields  of  application  has  been  the  main  thrust  of  this  research.  Previous 
two  reports  have  presented  the  mathematical  and  computational  developments  for  BVPs 
and  IVPs  in  hpk  framework  in  which  the  order  of  approximation  space  k  defining  global 
differentiability  of  order  (k  —  1)  has  been  pointed  out  as  an  independent  computational 
parameter  in  addition  to  h  and  p.  Successful  applications  of  this  new  hpk  framework  have 
been  presented  in  various  areas  of  continuum  mechanics  to  demonstrate  the  benefits  of 
using  this  framework  as  apposed  to  h,p  framework  used  currently  in  the  finite  element 
processes. 

The  thrust  of  the  work  completed  during  the  third  year  of  these  grants  have  been  in 
the  following  areas: 

(I)  Preliminary  research  towards  development  of  mathematical  models  for  fluid-solid 
interaction  and  associated  computational  infrastructure  in  hpk  mathematical  frame¬ 
work 

(II)  Preliminary  research  work  towards  development  of  concepts  for  a  priori  and  a  pos¬ 
teriori  error  estimation  in  h,  p,  k  mathematical  and  computational  framework 

(III)  Rediscretization,  moving  meshes  and  solution  mapping  strategies  and  associated 
computational  infrastructure  for  BVPs  and  IPVs  in  hpk  framework 

While  the  works  in  (I)  and  (II)  are  aimed  towards  preliminary  concepts  and  methodol¬ 
ogy  developments,  the  research  in  (III)  is  at  a  higher  maturity  level  and  has  already  been 
applied  to  many  BVPs  and  IVPs,  but  some  work  is  needed  to  streamline  the  algorithms 
for  applications  to  practical  problems  of  interest.  Details  for  each  of  the  three  areas  of 
research  are  presented  in  the  following. 

2  Fluid-solid  interaction 

The  preliminary  work  dealing  with  general  methodology  and  development  of  mathematical 
models  for  fluid-solid  interaction  problems  (BVPs  and  IVPs)  is  presented  in  this  section. 

2.1  Introduction 

The  thrust  of  this  research  is  to  develop  basic  concepts  and  methodologies  towards  de¬ 
velopment  of  a  new  mathematical  and  finite  element  computational  framework  based  on 
h,  p,  k  for  fluid-solid  interaction  BVPs  and  IVPs.  For  all  fluid-solid  interaction  problems, 
a  single  mathematical  model  describing  all  phases  of  the  physics  (i.e.  solid,  liquid,  or 
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gas)  is  essential  in  the  development  of  a  rigorous  mathematical  framework  that  natu¬ 
rally  translates  into  a  sound  computational  infrastructure  in  which  no  additional  ad-hoc 
adjustments  and  treatments  are  needed.  The  interface  behaviors  are  intrinsic  in  the  math¬ 
ematical  models  hence,  details  at  solid-liquid,  liquid-gas  or  solid-gas  interfaces  are  resolved 
naturally  and  accurately.  This  methodology  eliminates  the  need  for  constraint  equations 
at  the  interfaces  to  couple  the  behaviors  of  diverse  media  in  complex  interaction  prob¬ 
lems.  It  also  eliminates  inter-element  flux  problems  and  problem  dependent  treatments 
associated  with  moving  meshes  and  re-discretizations. 

While  the  mathematical  models  are  generally  derived  using  conservation  laws,  consti¬ 
tutive  equations,  and  equations  of  state,  the  specific  forms  of  the  governing  differential 
equations  (GDEs)  vary  significantly.  First-order  systems  derived  using  auxiliary  variables 
and  auxiliary  equations,  decomposition  of  total  stress  into  elastic  and  viscous  parts  for 
polymer  flows,  use  of  vorticity  and  stream  functions  are  a  few  examples  of  commonly 
used  approaches.  When  using  the  Galerkin  method  with  weak  form,  some  forms  of  GDEs 
may  be  advantageous  over  the  others.  The  current  mathematical  framework  based  on  h,  p 
finite  element  processes  employs  local  approximations  of  class  C°  in  space  as  well  as  in 
time.  Convergence  of  the  numerical  errors  is  measured  in  terms  of  Lebesgue  measures  or 
■^2-type  global  norms.  Consequently,  the  differentiability  requirements  may  be  violated 
at  sets  of  measure  zero  in  the  discretizations.  This  results  in  flexibility  but  at  the  expense 
of  local  (inter-element)  accuracy  and  local  convergence.  The  proposed  framework  utilizes 
the  theory  of  continuous  and  differentiable  functions  in  conjunction  with  Sobolev  spaces 
and  hence  it  is  possible  to  incorporate  the  desired  global  differentiability  in  the  design  of 
the  computational  process  due  to  k,  the  order  of  the  approximation  space. 

In  this  research  it  is  demonstrated  that  the  GDEs  in  their  strong  forms  (higher  order 
systems)  are  highly  meritorious  over  all  others.  Additionally,  use  of  fundamental  variables 
such  as  density,  velocities,  temperature  is  preferred  over  others  (mathematical  quantities) 
in  that  such  variables  naturally  lead  to  simplicity  in  defining  boundary  and  initial  condi¬ 
tions.  Use  of  the  strong  form  of  GDEs  also  eliminates  redundancies  and  inconsistencies 
that  may  arise  when  constructing  local  approximations.  Higher  order  systems  of  GDEs, 
however,  require  approximations  in  higher-order  spaces  that  have  many  additional  benefits 
[1-4]  as  shown  subsequently. 

In  Section  2.2,  we  present  a  brief  summary  of  currently  used  finite  element  approaches, 
followed  by  the  development  of  mathematical  models,  mathematical  and  computational 
approach  in  Sections  2. 3-2. 5. 

2.2  Background 

2.2.1  Methods  of  Approximation 

In  designing  finite  element  processes  one  utilizes  mathematical  models  to  construct  inte¬ 
gral  forms  [5-7].  For  example,  consider  the  operator  equation  A<p  —  f  —  0  in  D,  whore 
A  is  a  differential  operator,  <p  is  the  dependent  unknown,  and  /  is  the  data.  We  assume 
that  the  solution  ip  exists  and  is  unique.  We  construct  /n(A</?  —  f)v  dfl  =  0.  Here  v  is 
an  arbitrary  weight  (test)  function  that  vanishes  on  the  part  where  =  <p0-  Utilizing 
this  integral  form  we  can  distinguish  between  the  following  weighted-residual  methods: 
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1.  Petrov-  Galerkin  method  (PGM):  In  this  method  the  weight  function  is  an  indepen¬ 
dent  function  v  =  ip  ^  5ifh,  with  v  =  0  on  Ti  if  p  =  <po  on  I\. 

2.  Galerkin  method  ( GM):  If  <p/,  is  an  approximation  of  <p,  the  Galerkin  method  seeks 
<P/t  such  that  Jn(A(ph  —  f)v  dQ  =  0  with  v  =  5(ph-  This  can  be  written  as  Bg(v,  iph)  — 
lg{v). 

3.  Galerkin  method  with  weak  form  (GMWF):  In  this  approach,  we  set  v  =  5iph 
and  transfer  some  of  the  differentiation  from  <p/,  to  5iph  to  obtain  Bgw(6<ph,  iph)  = 

l gw  * 

4.  Least  squares  method  (LSM):  In  this  approach  we  define  the  residual  E  by  E  = 
A<ph  —  /,  construct  a  functional  I(<Ph)  —  (E,  E) q,  and  determine  iph  by  minimiz¬ 
ing  I(tph)-  Thus,  SI(iph)  —  0  provides  the  necessary  condition  from  which  iph  is 
determined,  and  52I{^h)  >  0  provides  the  sufficient  condition. 

2.2.2  Variationally  Consistent  Finite  Element  Processes 

The  methods  of  approximation  described  above  are  classical  methods  in  which  the  domain 
of  definition  Q  is  not  discretized.  The  major  drawback  of  these  methods  is  that  it  is  very 
difficult  to  construct  global  approximation  functions  for  arbitrary  domains  and  boundary 
conditions  [7-10].  In  FEM,  wc  utilize  these  methods  over  a  discretized  domain  LtT  =  Uefie 
in  which  Qe  is  a  closure  of  Lle,  an  element.  In  this  approach,  the  global  approximation 
<p/i  =  U eipeh)  where  <p\  is  the  restriction  of  <p/,.  over  f!e,  and  is  referred  to  as  the  local 
approximation.  While  the  local  differentiability  of  q>h  can  be  controlled  by  the  degree  of 
local  approximation  p,  the  global  differentiability  of  <ph  depends  entirely  on  the  order  k  of 
the  approximation  space.  Thus,  k  is  an  independent  parameter  in  finite  element  processes 
in  addition  to  h  and  p.  The  need  for  global  differentiability  arises  from  the  requirements 
of  continuity  and  differentiability  of  <p  over  which  is  intrinsic  in  the  physics  (see  Figure 
1  and  Sections  3.2  and  3.3).  The  issues  of  non-smooth  data  and  applied  disturbance  often 
leading  to  singular  BVP  and  IVP  can  also  be  easily  resolved  in  h,p ,  k  framework  [11], 

In  the  GM,  GMWF  and  PGM,  we  convert  a  system  of  differential  equations  into  an 
algebraic  system  using  integral  forms  and  the  approximation  < ph-  In  the  current  math¬ 
ematical  framework,  we  must  establish,  on  a  problem  by  problem  basis,  as  to  which 
possible  choices  of  h  and  p  ensure  well-posedness.  This  issue  is  resolved  in  the  current 
research  by  establishing  a  link  between  the  integral  forms  and  the  elements  of  the  calcu¬ 
lus  of  variations.  Calculus  of  variations  deals  with  extrema  of  functionals  which  are  often 
integrals. 

Let  I(iph)  be  a  functional.  Then  based  on  calculus  of  variations  [5-7],  51(<ph)  =  0 
provides  the  necessary  condition  from  which  < ph  is  determined  and  52/(<p/l)  provides  an 
extremum  principle  or  sufficient  condition  that  establishes  uniqueness  of  tph  obtained 
from  81{tph)  =  0.  One  can  show  that  a  <p/,  that  yields  a  unique  extrema  of  7(<p/()  is  also  a 
unique  solution  of  the  differential  equation  Aq>h  —  f  =  0,  which  is  often  referred  to  as  the 
Euler  equation.  Thus,  we  have  a  link  between  the  solution  iph  of  the  differential  equation 
and  calculus  of  variations.  We  note  that  integral  forms  in  the  GM,  GMWF,  and  PGM 
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are  equivalent  to  Sl(<Ph)  =  0.  In  these  integral  forms,  existence  of  I (if/,)  is  not  always 
guaranteed.  Even  when  I(<ph)  exists,  J2/(yp^)  may  not  yield  a  unique  extremum  principle. 

To  distinguish  well-posedness  of  the  integral  forms,  the  concepts  of  variational  consis¬ 
tency  (VC)  and  variational  inconsistency  (VIC)  are  introduced.  If  there  exists  a  functional 
I(<fh)  such  that  5I((fih)  =  0  and  62I(iph)  yields  a  unique  extremum  principle  then  the  inte¬ 
gral  form  under  consideration  is  termed  variationally  consistent.  Obviously,  VC  integral 
forms  always  ensure  unique  iph-  Integral  forms  that  are  not  VC  are  called  variationally 
inconsistent.  In  such  cases  the  resulting  computations  may  not  be  unconditionally  stable 
or  may  fail  all  together.  It  has  been  shown  by  Surana  et.  al  [1-4]  that  the  concepts 
of  VC  and  VIC  permit  us  to  decide  whether  integral  forms  would  yield  unconditionally 
positive-definite  algebraic  systems  and  hence  unconditional  numerical  stability.  This  is  an 
issue  of  paramount  importance  in  adaptive  processes  in  which  one  often  begins  with  very 
coarse  discretizations  and  very  low  p-levels.  Variationally  inconsistent  integral  forms  in 
such  cases  may  not  even  permit  computations  without  the  use  of  ad-hoc  approaches  and 
upwinding  methods.  The  PGM  and  GM  are  always  VIC  regardless  of  the  differential  op¬ 
erators;  GMWF  is  variationally  consistent  when  the  bilinear  functional  is  continuous  and 
satisfies  the  so-called  Inf-Sup  condition  which  is  only  possible  for  self-adjoint  operators; 
LSM,  on  the  other  hand,  is  always  VC  for  all  three  types  of  differential  operators.  These 
conclusions  are  significant  and  helpful  in  assessing  the  current  finite  element  computa¬ 
tional  technology  and  associated  deficiencies  as  well  as  in  designing  a  new  computational 
framework  in  which  the  computational  processes  always  remain  numerically  stable  and 
nondegenerate. 

2.2.3  Mathematical  Models  and  Currently  Used  Finite  Element  Methods 

The  mathematical  models  for  fluid-solid  interaction  problems  are  diverse.  The  mathe¬ 
matics  of  performing  finite  element  computations  depends  on  whether  the  models  are 
continuum  based.  In  the  proposed  research  we  only  consider  fluid-solid  interaction  math¬ 
ematical  models  that  are  based  on  continuum  mechanics.  The  solid  medium  could  be 
elastic,  elasto-plastic,  visco-elastic  or  even  energetic  material.  The  fluid  can  be  an  incom¬ 
pressible  Newtonian,  Generalized  Newtonian  or  a  polymeric  liquid  with  elasticity  as  well 
as  viscosity.  The  fluid  medium  can  also  be  a  gas  with  complex  properties  requiring  real 
gas  models  and  temperature  dependent  transport  properties.  The  combined  fluid-solid 
problems  may  be  non-isothermal  with  large  motions,  large  strain  rates,  moving  interfaces, 
and/or  moving  boundaries.  In  addition,  the  problems  may  involve  chemical  kinetics  if 
the  materials  are  energetic  or  reacting.  Towards  obtaining  numerical  solutions  of  such 
interaction  problems,  we  present  a  comprehensive  review  of  currently  used  mathematical 
and  finite  element  models.  This  is  followed  by  details  of  the  concepts  developed  in  the 
present  methodology. 

The  mathematical  models  of  fluid-solid  interactions  often  result  in  non-linear  partial 
differential  equations  in  space  coordinates  (BVP)  or  space  cordinates  and  time  (IVP).The 
most  prevalent  finite  element  computational  methodologies  for  obtaining  their  numer¬ 
ical  solutions  is  based  on  GMWF  which  results  in  integral  forms  that  are  VIC.  Such 
integral  forms  yield  computational  processes  that  are  not  unconditionally  stable  or  non¬ 
degenerate.  Thus,  one  must  establish  for  each  problem,  the  ranges  of  computational  and 
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physical  parameters  for  which  the  computations  will  remain  stable.  This  often  leads  to 
unreasonable  and  impractical  dicretizations  and  p-levels.  The  desire  to  compute  with 
coarser  dicretizations  and  lower  p-levels  has  lead  to  the  developments  of  upwinding  meth¬ 
ods  [12-26].  Amongst  many  other  drawbacks,  the  most  disturbing  aspect  of  such  methods 
is  that  they  are  problem  dependent.  Hence,  they  cannot  be  considered  as  a  general  com¬ 
putational  methodology  for  all  BVP  and  IVP.  In  the  proposed  research  we  only  consider 
VC  integral  forms.  We  therefore,  eliminate  the  need  and  use  of  problem  dependent  ap¬ 
proaches  all  together. 

The  fluid-solid  interaction  problems  most  naturally  can  be  viewed  as  initial- value  prob¬ 
lems  (IVP).  Current  computational  approaches  for  IVP  can  be  broadly  classified  into  two 
groups:  space-time  decoupled  methods  and  space-time  coupled  methods.  Space  time 
methods  with  time  marching  (as  opposed  to  space-time  meshes)  offer  many  advantages 
[4,  27].  Most  significant  features  of  space-time  coupled  methods  are:  (i)  concurrent  de¬ 
pendence  of  the  solution  on  space  and  time,  (ii)  possibility  of  making  mathematical  classi¬ 
fications  of  space-time  differential  operators  and  hence  developing  a  general  mathematical 
framework  for  all  IVP,  (iii)  possibility  of  assessing  VC  or  VIC  of  space  time  integral  forms, 
and  (iv)  VC  integral  forms  lead  to  unconditionally  non-degenerate  computations  (no  is¬ 
sues  of  stability).  None  of  these  features  exist  in  space-time  decoupled  methods.  In  the 
present  work,  we  only  utilize  space-time  integral  forms  that  are  space-time  variationally 
consistent. 

The  mathematical  models  of  fluid-solid  interaction  problems  can  be  established  either 
using  Lagrangian  approach  or  Eulerian  approach.  In  the  Lagrangian  approach,  most  com¬ 
monly  used  in  solid  mechanics,  all  quantities  of  interest  are  measured  in  a  fixed  frame  of 
reference  embedded  in  the  body.  In  contrast,  in  the  Eulerian  description  all  quantities  of 
interest  are  measured  and  referred  to  a  reference  frame  fixed  in  space  and  material  parti¬ 
cles  move  through  the  space.  In  the  currently  used  approaches  for  fluid-solid  interaction 
problems,  each  medium  of  the  process  is  mathematically  modeled  and  discretized  using 
the  most  suitable  strategy  for  each  medium  and  then  interfaced  through  constraint  equa¬ 
tions  to  the  neighboring  media.  In  this  approach:  (1)  The  interface  behavior,  which  is 
often  of  great  interest  is  highly  dependent  on  the  nature  of  constraint  equations  describing 
the  coupling.  Their  legitimacy  in  relation  to  the  true  physical  behavior  may  be  difficult 
to  ascertain.  (2)  Often  the  variables  from  each  medium  at  the  interface  boundaries  are 
different.  This  results  in  significant  difficulties  in  describing  constraint  equations  at  the 
interfaces.  (3)  The  most  disturbing  feature  of  this  approach  is  that  the  integration  of 
the  different  media  of  the  process  at  the  interfaces  is  only  established  numerically  after 
discretizations.  (4)  Inaccuracies  in  the  constraint  equations  at  the  interfaces  may  also  con¬ 
taminate  the  solutions  in  the  entire  media.  Such  approaches  are  problem  dependent  and 
can  not  viewed  as  a  general  methodology  for  all  fluid-solid  interaction  problems.  In  the 
approach  presented  here,  all  media  of  fluid-solid  interaction  problems  are  mathematically 
modeled  under  one  mathematical  umbrella.  This  approach  completely  eliminates  the  ne¬ 
cessity  of  applying  interface  constraint  equations  but  requires  the  development  of  rigorous 
mathematical  models.  Details  of  this  approach  are  described  in  Section  2.3.  In  summary, 
the  mathematical  models  in  Lagrangian  frame  of  reference  for  solids  and  those  for  fluid 
in  Eulerian  frame  of  reference  with  constraint  equations  at  the  interfaces  describing  their 
coupling  remains  the  current  methodology  for  fluid-solid  interaction  processes. 
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2.2.4  Moving  Fronts,  Boundaries,  Meshes,  and  Rediscretizations 

Accurate  and  practical  numerical  simulations  of  fluid-solid  interaction  problems  may  re¬ 
quire  a  representation  of  moving  interfaces,  fronts  and/or  moving  boundaries  that  can 
be  tracked  accurately.  The  discretizations  in  other  regions  are  kept  intact  but  moved  to 
different  locations.  This  requires  a  mathematical  infrastructure  in  which  front  tracking 
and  moving  mesh  capability  is  transparent  during  computations.  In  other  instances  a 
total  re-discretization  may  be  required.  The  mathematical  framework  to  support  these 
features  are  of  critical  significance  in  interaction  problems.  It  requires  the  results  of  one 
discretization  to  be  mapped  onto  a  new  dicretization.  It  is  well  known  that  in  the  h,p 
framework  based  on  global  approximations  of  class  C°,  the  mapping  of  solutions  from  one 
discretization  to  another  may  suffer  severe  damage  in  the  solution  derivatives,  referred 
to  as  the  flux  problem  in  the  published  literature  (discontinuity  of  normal  derivative  to 
the  inter-element  boundaries,  see  Figure  1).  Many  problem  dependent  remedies  are  avail¬ 
able  [27]  to  circumvent  these  difficulties  to  some  extent.  Due  to  the  C°  nature  of  the 
global  approximation  used,  the  inter-element  flux  discontinuity  remains.  Thus,  within 
h,p  framework  with  C°  global  approximations,  a  general  methodology  for  moving  fronts, 
moving  meshes  and  re-discretization  that  is  free  of  flux  problem  is  not  possible. 

Moving  interfaces,  fronts  and  boundaries  can  be  tracked  using  the  level  set  theory.  The 
h,  p,  k  framework  maintains  sharpness  of  the  fronts  during  all  stages  of  the  evolution.  Since 
the  global  differentiability  of  the  approximation  is  controlled  through  k,  flux  problems  are 
totally  absent  in  the  present  approach.  Mapping  of  the  solutions  from  one  discretization  to 
another  with  desired  global  differentiability  and  measures  of  their  accuracy  remain  totally 
transparent  (See  Section  4  on  Rediscretization, moving  meshes  and  solution  mapping). 

2.3  Development  of  Mathematical  Models 

This  is  one  of  the  most  important  phases  of  the  present  research.  We  recognize  that  the 
fluid  description  in  Eulerian  frame  of  reference  is  essential  due  to  large  motions  of  the 
fluid  particles.  In  the  present  work,  we  construct  the  mathematical  model  of  the  entire 
fluid-solid  system  in  the  Eulerian  frame  of  reference. 

The  mathematical  model  for  the  fluid  medium  are  constructed  using  the  Eulerian  frame 
of  reference  and  conservation  of  mass,  momentum  and  energy,  constitutive  equations, 
and  equations  of  state.  The  fluid  may  be  Newtonian,  generalized  Newtonian  (power  law, 
Carreau-Yasuda,  Bingham  etc.),  or  viscoelastic  liquid  (Maxwell,  Oldroyd,  Giesekus,  PTT, 
etc.),  or  a  gas  with  complex  constitutive  models  and  temperature  dependent  transport 
properties.  The  process  can  be  non- isothermal  with  chemical  kinetics.  The  variables  in  the 
formulation  may  include  density,  velocities,  pressure,  stresses,  and  temperature,  among 
others.  For  the  equations  of  state  describing  relationships  between  density,  pressure  and 
temperature,  we  consider  ideal  gas  as  well  as  real  gas  models  such  as  Van  der  waals, 
Redlich-Kwong,  Beattie-Bridgeman,  Benedict-Web-Rubin,  and  virial  models.  Transport 
properties  such  as  viscosity,  thermal  conductivity,  specific  heat  will  be  considered  to  be 
temperature  and  possibly  strain-rate  dependent.  Models  such  as  power  law,  Southerland 
law  are  utilized  for  this  purpose. 

For  a  linear  or  non-linear  elastic  solid  medium  as  well  as  the  solid  medium  undergo- 
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ing  large  deformation  and  strain  rates  with  plasticity,  the  mathematical  models  are  also 
derived  in  the  Eulerian  frame  of  reference.  While  the  basic  conservation  laws:  continu¬ 
ity,  momentum  and  energy  remain  the  same  as  in  the  case  of  fluids,  some  modifications 
are  needed  for  solids.  Additionally,  the  constitutive  equations  for  solids  require  more 
careful  considerations.  The  solid  medium  may  be  elastic,  elasto-plastic,  viscoelastic,  het¬ 
erogeneous  (i.e. , laminated  composite)  or  even  energetic.  The  Cauchy  stress  tensor  er  is 
decomposed  into  deviatoric  and  hydrostatic  parts.  With  this  decomposition,  deviatoric 
stress  components  in  solids  have  the  same  physical  meaning  as  viscous  stress  in  case  of 
fluids  after  using  the  Stokes  hypothesis.  Density,  velocities,  pressure,  deviatoric  part  of 
Cauchy  stress,  and  temperature  are  the  variables  of  choice  in  the  mathematical  models 
of  both  fluid  and  solid  media. 

Use  of  velocities  in  the  conservation  laws  requires  that  the  constitutive  equations  for 
the  solid  medium  be  derived  in  terms  of  strain  rates  as  opposed  to  strains  used  in  the 
Lagrangian  frame  of  reference.  For  this  purpose,  we  employ  rate  constitutive  models 
based  on  stress  tensors  of  Jaumann,  Truesdell,  Green-Nagdi,  and  convected  stress  tensor 
[28-31]: 


i  Do 
~Dt 


=  jC  :  D; 


tDo 

Dt 


=  ‘C  :  D; 


gnDo 

Dt 


=  gnC  :  D; 


cDo 

Dt 


CC  :  D 


where  with  superscripts  j,  t,  gn  and  c  are  Jaumann,  Truesdell,  Green-Nagdi,  and 
convected  rates  of  Cauchy  stress  a.  Here  JC,  tC,  9nC  and  CC  are  the  corresponding 
material  tangent  moduli  tensors  and  D  is  the  symmetric  part  of  the  rate  of  deformation 
tensor.  The  rate  constitutive  equations  are  invariant  of  rotations.  Note  that  the  choice 
of  one  rate  equation  over  other  is  material  dependent  and  there  is  no  consensus  on  their 
choice.  In  the  present  research,  we  consider  all  four  rate  equations  and  investigate  their 
ranges  of  applications  for  various  materials.  When  the  solid  medium  is  compressible  (as 
in  case  of  high  strain  rate  plasticity)  appropriate  equations  of  state  (i.e.,  [32]  and  others) 
arc  utilized  to  define  relationships  between  density,  pressure  and  temperature. 

When  these  mathematical  models  are  discretized  and  solved  numerically,  the  inter¬ 
action  of  the  solid  medium  with  the  fluid  medium  at  the  interfaces  is  intrinsic  in  the 
mathematical  models  and  hence  completely  eliminates  the  use  of  constraint  equations  at 
the  interfaces.  Since  the  mathematical  models  are  in  Eulerian  frame  of  reference,  the  com¬ 
putations  for  these  models  are  on  a  fixed  discretization  hence  eliminating  re-discretizations 
and  mapping  of  solutions  as  integral  part  of  the  numerical  computational  methodology. 


2.4  Higher-Order  Global  Differentiability  of  Approximation  Func¬ 
tions 

Theoretical  solutions  of  the  mathematical  models  may  be  of  higher-order  global  differen¬ 
tiability  in  space  as  well  as  in  space  and  time.  When  considering  numerical  solutions  of 
these  governing  differential  equations  using  FEM,  the  higher-order  global  differentiable 
approximations  in  space  as  well  as  in  space  and  time  are  meritorious  over  solutions  of 
class  C°  [1-4]. 

Surana  et.  al  [1-4]  have  shown  that  the  order  k  of  the  approximation  space  Hk’p 
establishes  global  differentiability  of  order  (A;  —  1)  of  the  numerical  solutions,  and  is  an 
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independent  parameter  in  all  finite  element  computations  in  addition  to  h  and  p.  They 
have  developed  a  k-ve rsion  of  the  finite  element  method  and  associated  hp,  pk  and  hpk 
processes.  Authors  also  shown  that  even  the  most  prudent  hp-adaptive  processes  cannot 
alter  the  order  k  of  the  approximations.  Adaptive  control  of  higher-order  global  differ¬ 
entiability  allows  control  of  continuity  and  differentiability  of  any  quantity  of  interest  in 
the  design  of  a  computational  process.  The  higher-order  globally  differentiable  approxi¬ 
mations  in  space  and  time  for  initial- value  problems  (I VP)  require  the  use  of  Hk,p  spaces 
in  which  k  =  (Aq,  k2)  and  p  =  (pi,p2)  with  pi  >  2 Aq  -  1  and  p2  >  2k2  —  1;  here  Aq  and 
k2  denote  the  orders  of  the  approximation  spaces  in  spatial  coordinates  and  time,  respec¬ 
tively,  and  p\,p2  are  the  corresponding  degrees  of  approximations.  Higher-order  global 
differentiability  of  the  approximations  is  necessary  when  the  theoretical  solution  is  of  a 
higher  class  and  therefore  has  higher-order  global  derivatives.  In  the  case  of  singular  data, 
boundary  or  initial  conditions  and  their  approximations  in  Hk,v  space  using  interpolants 
permit  one  to  resolve  the  desired  behaviors  as  a  suitable  k  is  adapted  [11],  Use  of  Hk’p 
scalar  product  spaces  for  space-time  approximations  has  been  extensively  investigated  by 
Surana  et.  al  [1-4,  11]  and  they  have  shown  that  very  high  accuracy  is  achievable. 

For  illustrative  purposes,  we  consider  the  steady  state  convection-diffusion  equation 
(BVP)  for  Peclet  number  of  100  shown  in  Figure  1.  The  theoretical  solution  is  of  class 
C°°.  Numerical  solutions  of  classes  C-7,  j  =  0,..,4  (corresponding  to  k  =  1,...,5)  are 
computed  using  least  squares  finite  element  formulation  (VC)  for  a  uniform  discretization 
of  ten  p- version  ID  elements.  In  most  engineering  applications,  at  least  the  first  derivative 
of  the  solution  is  of  considerable  interest.  Figure  1  shows  diph/dx  versus  x  at  and  in  the 
neighborhood  of  x  =  0.9  (inter-element  boundary)  for  the  solutions  of  various  classes  and 
a  comparison  with  the  theoretical  solution.  Discontinuity  of  dph/dx  for  the  solutions  of 
class  C°  (k  =  1)  at  x  =  0.9  and  its  influence,  causing  oscillations  in  the  neighboring 
elements  is  clearly  observed.  The  discontinuity  of  dphjdx  is  also  observed  at  all  other 
inter-element  boundaries,  however  the  magnitude  decreases  moving  from  x  =  1.0  toward 
x  =  0.0  due  to  decreasing  magnitude  of  dpjdx  (theoretical  value).  Solutions  of  classes 
C 1  to  C4  show  dramatically  improved  behavior  of  dphjdx  with  lower  degrees  of  freedom 
(dofs)  compared  to  those  used  to  obtain  the  C°  solution.  Solution  of  class  C4  is  almost 
in  perfect  agreement  with  the  theoretical  solution  for  only  108  dofs  compared  to  120  dofs 
used  to  obtain  the  C°  solution  that  is  grossly  in  error.  Progressive  reductions  in  total 
dofs  for  progressively  higher  classes  of  approximations  are  more  dramatic  for  2D  and  3D 
problems.  We  note  that  the  L2  norm  of  the  error  in  the  first  derivative  for  C4  solution  is 
two  orders  of  magnitude  lower  than  that  for  C°  solution.  Ability  of  the  approximations 
of  higher  classes  in  approaching  theoretical  solutions  efficiently  is  clearly  observed.  In 
this  research,  we  propose  to  construct  the  space  and  space-time  approximations  of  higher 
order  global  differentiability  in  the  Hk  p  spaces  for  obtaining  finite  element  solutions  of 
the  mathematical  models. 

2.5  h,p,k  Finite  Element  Processes:  VC  Integral  Forms 

The  governing  equations  of  non-linear  processes  in  the  fluid-solid  interaction  considered 
here  are  naturally  non-linear  partial  differential  equations.  In  case  of  IVP  all  quantities  of 
interest  are  simultaneously  dependent  on  space  coordinates  as  well  as  time.  We  solve  these 
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GDEs  for  both  BVP  and  IVP  using  h,  p,  k  finite  element  processes  in  which  approximations 
in  space  as  well  as  in  space  and  time  are  of  higher-order  global  differentiability.  The  space 
as  well  as  space-time  integral  forms  are  always  variationally  consistent.  The  significant 
features  of  this  methodology  are  summarized  below. 

•  For  BVPs,  finite  element  processes  based  on  variationally  consistent  integral  forms 
are  considered  so  that  unconditional  stability  of  the  resulting  computational  pro¬ 
cesses  is  ensured. 

•  For  IVPs,  space-time  coupled  FE  processes  [4]  are  used  as  opposed  to  space-time 
decoupled  approaches  [27].  Space-time  coupled  methods  are  consistent  with  the 
physics  of  the  IVPs. 

•  Only  for  space-time  coupled  finite  element  processes  with  space-time  integral  forms 
it  is  possible  to  establish  a  link  between  the  integral  forms  and  the  calculus  of 
variations.  Surana  et.  al  [4]  have  shown  that  space-time  variationally  inconsistent 
integral  forms  yield  computational  processes  in  which  the  coefficient  matrices  are 
always  non-symmetric.  Such  computational  processes  are  only  conditionally  stable 
and  may  even  become  totally  degenerate  for  some  choices  of  h ,  p  and  k.  Space-time 
variationally  consistent  integral  forms  always  yield  symmetric  coefficient  matrices 
that  are  positive-definite  and  the  computational  models  are  unconditionally  stable 
and  non-degenerate  regardless  of  the  choice  of  h ,  p,  k  as  well  as  the  parameters  of  the 
physics.  The  space-time  differential  operators  in  the  mathematical  models  of  fluid- 
solid  interactions  are  naturally  non-linear.  Surana  et.  al  [4]  have  shown  that  for 
such  differential  operators  only  the  space-time  integral  forms  derived  using  space- 
time  least  squares  processes  are  VC  and  hence  are  computationally  meritorious  over 
all  others.  Therefore,  these  are  be  utilized  in  the  present  research. 

•  The  space-time  least  squares  processes  are  also  completely  free  of  numerical  disper¬ 
sion  [33]  for  appropriate  choices  of  h,p  and  k.  This  feature  permits  computations 
that  are  time-accurate  for  IVPs. 

•  Commonly  encountered  limitations  of  CFL  number  and  stability  issues  in  space-time 
decoupled  methods  are  totally  absent  in  the  present  methodology  due  to  variational 
consistency  of  the  integral  forms. 

•  When  GMWF  is  used  to  obtain  numerical  solutions,  the  resulting  integral  forms 
are  VIC,  and  it  becomes  necessary  to  use  upwinding  techniques  such  as  SUPG, 
SUPG/DC,  SUPG/DC/LS  [12-26]  and  so  on  to  stabilize  the  computations.  In  the 
current  methodology,  VC  of  the  least  squares  integral  forms  eliminates  the  use  of 
problem  dependent  upwinding  methods  all  together.  This  is  a  significantly  advan¬ 
tageous  characteristic  of  the  computational  methodology  presented  here. 

•  Surana  et.  al  [4]  and  others  [27]  have  shown  that  space-time  marching  techniques 
are  advantageous  over  the  space-time  meshes  in  terms  of  computational  efficiency 
as  well  as  control  over  the  solution  errors  during  evolution  of  the  solution.  The 
numerical  solution  is  computed  for  one  strip  at  a  time,  and  only  when  the  desired 
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solution  accuracy  is  achieved  then  marched  in  time  until  the  target  time  is  reached. 
In  space-time  meshes,  the  problem  size  is  significantly  larger  than  that  in  the  space- 
time  strip  approach,  and  since  the  solution  is  computed  simultaneously  for  all  values 
of  time,  error  control  as  in  space-time  strip  is  not  possible. 

•  While  much  is  published  in  the  numerical  analysis  of  fluid-solid  interaction  problems, 
the  solutions  independent  of  the  computational  parameters  (referred  to  as  mesh 
independent  solutions)  are  generally  not  available.  In  h,  p,  k  framework  proposed 
here,  numerical  simulations  are  sought  that  are  independent  of  h  and  p  for  a  chosen  k. 
The  proximity  of  such  solutions  to  the  theoretical  solutions  are  established  through 
appropriate  measures. 

•  For  a  numerical  solution  of  a  BVP  to  be  valid  it  must  correspond  to  the  stationary 
state  of  the  corresponding  IVP.  In  many  instances  in  polymer  flows  in  which  the 
constitutive  equations  are  highly  complicated  and  are  empirical  (Maxwell,  Oldroyd, 
Giesekus,  PTT,  models  etc.).  This  aspect  is  essential  in  establishing  the  legitimacy 
of  the  solutions  of  the  BVP. 

•  Time  accurate  evolutions  for  fluid-solid  interaction  problems  are  generally  not  re¬ 
ported.  Numerical  solutions  of  such  processes  available  in  the  literature  are  obtained 
using  GMWF  (VIC)  with  upwinding  techniques.  In  the  present  work,  we  establish 
time  accurate  evolutions  of  fluid-solid  interaction  problems. 

•  In  the  following  we  present  the  basic  steps  of  LSM  for  non-linear  BVP  and  IVP  [3-4] 
that  are  encountered  in  fluid-solid  interactions.  Let  Aip  —  /  =  0  be  a  BVP  or  an 
IVP  in  n,  tfh  be  an  approximation  of  p  and  Aph  —  f  =  E  be  residual  in  0.  Then 
we  construct  a  functional, 

I(V„)  =  (E,E)  (1) 

and  the  necessary  condition  for  its  minimum  is 

6I(<ph)  =  (E,6E)  =  g(<ph)  =  0  (2) 

The  second  variation  yields 

82I(ph)  =  (8E,8E)  +  (E,82E)  (3) 

We  note  that  52I(iph)  is  not  unconditionally  greater  than  zero  and  does  not  represent 
a  unique  extremum  principle  unless  (E,82E)  in  (3)  is  neglected  [3-4,33],  which 
amounts  to  correcting  the  search  direction.  However,  we  still  seek  the  solution  of 
(2).  Thus  we  use 

82I(ph)  =  (5E,8E)  (4) 

Note  that  82I{<ph)  defined  by  (4)  is  always  greater  than  zero;  hence,  a  unique  ex¬ 
tremum  principle  and  a  <ph  obtained  from  (2)  minimizes  I(<ph)  in  (1).  The  system 
of  non-linear  equations  in  (2)  are  solved  using  Newton’s  method  with  line  search. 
Let  iph°  be  a  starting  solution,  then,  we  have  the  following: 

<Ph  =  <Ph°  +  aA(ph0-,  A  (ph  =  -\82I{vh)]~^h{g{<-Ph)}vi  (5) 
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The  factor  a  is  determined  such  that  I(<Ph)  <  / (v3/!)-  We  remark  that  the  LSM  based 
on  (1),  (2),  (3),  and  (5)  is  VIC,  while  the  LSM  based  on  (1),  (2),  (4),  and  (5)  is  VC. 
Rationale  of  deriving  (4)  from  (3)  has  been  established  in  [3,34],  We  consider  the 
steady  state  Burgers  equation  in  (0,1)  with  y>(0)=1.0  and  </?(l)=0.0.  Figure  2  shows 
numerical  solutions  obtained  using  VC  and  VIC  least  squares  formulations.  At  p  = 
7,  VIC  formulation  produces  spurious  solution,  though  criterion  for  the  convergence 
of  the  Newton’s  method  with  line  search  is  satisfied  and  the  least  squares  functional 
is  of  the  same  order  of  magnitude  as  in  case  of  VC  integral  form.  Additionally, 
convergence  of  the  iterative  solution  method  slows  down  significantly  in  case  of  a  VIC 
integral  form  for  all  p  levels,  requiring  over  50  iterations,  whereas  VC  formulation 
always  converges  in  less  than  10  iterations  for  all  p  levels.  This  example  demonstrates 
the  potential  of  failure  of  VIC  integral  form  for  some  choice(s)  of  computational 
parameters  (p  level  in  this  case).  Consequences  of  VIC  integral  forms  in  GAL/WF 
are  far  more  serious  if  upwinding  methods  are  not  used. 

2.6  Preliminary  Results 

Figure  3  shows  preliminary  results  of  evolutions  for  ID  axial  stress  wave  propagation  in 
media  containing  bimaterial  interfaces  obtained  using  computational  methodology  dis¬ 
cussed  in  this  research.  The  mathematical  models  for  both  fluid  and  solid  media  are  in 
Eulerian  frame  of  reference.  Interface  coupling  equations  and  upwinding  methods  are 
not  used  at  all.  The  integral  forms  are  space-time  variationally  consistent,  hence  ensure 
unconditionally  stable  computations  during  all  stages  of  the  evolution.  In  all  cases  wave 
propagation,  transmission,  reflection  and  interface  behavior  is  simulated  easily.  Interface 
behaviors  are  oscillation  free  and  wave  motions  are  in  conformity  with  the  impedences  of 
the  media  and  their  transport  properties. 
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Figure  1  dtp/dx  versus  x  for  steady  state  convection  diffusion  equation  : 

Solutions  of  various  classes  (  exploded  view  at  x  =  0.9,  an  inter-element  boundary  ) 


0.6 


Solution  for  VC  integral  form  at  p  =  7  |ii 


- . . - . . . 

. \ 
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Spurious  solution  for  VIC  integral  form  at  p  =  7 
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Figure  2  Solution  <p  versus  x  for  steady  state  Burgers  equation  for  different  /7-levels  : 
VC  and  VIC  formulations 
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Wave  Propagation  :  Time  Evolutions 


Medium  1  :  Ml 


Medium  2  :  M2 


■> 


Ml  :  Newtonian  fluid 
M2  :  Soft  solid 
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Ml  :  Newtonian  fluid 
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Ml  :  Dense  polymer 
M2  :  Steel 


Ml  :  Elastic  solid  1 
M2  :  Elastic  solid  2 


Figure  3  Stress  wave  propagation,  transmission  and  reflection  in  bimaterial  media. 
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3  A  priori  and  a  posteriori  error  estimation  :  Pre¬ 
liminary  work 

3.1  Introduction 

In  recent  works,  Surana  et.al.  [1-4]  have  shown  that  order  k  of  the  approximation  space 
defines  global  differentiability  of  order  (k  —  1)  of  the  approximations  and  is  an  indepen¬ 
dent  parameter  in  all  finite  element  computations.  Authors  have  established  k-  version 
of  finite  element  method  and  associated  hk,pk ,  hpk  processes  in  addition  to  h,  p- versions 
and  h,p,hp  processes  used  currently.  Under  the  AFOSR  grants  F  49620-03-1-0298  and 
F  49620-03-1-0201,  Surana  et.al.  [1-4]  proposed  and  developed  hpk  mathematical  and 
computational  framework  for  finite  element  processes  for  solving  BVP  and  IVP.  Due  to 
the  fact  that  h,  p  and  k  are  independent  parameters  in  the  new  mathematical  framework, 
all  quantities  of  interest  in  the  finite  element  computations  within  this  framework  are 
naturally  dependent  on  h,  p  and  k  as  opposed  to  only  h  and  p. 

A-priori  and  a  posteriori  error  estimates  are  intrinsically  important  aspects  of  a  math¬ 
ematical  framework  for  finite  element  computations.  A  priori  and  a  posteriori  error  esti¬ 
mates  derived  and  used  currently  are  strictly  based  on  h  and  p.  The  global  approximations 
are  almost  always  assumed  to  be  of  class  C°.  The  purpose  of  this  research  is  to  examine 
the  methodology  for  a  priori  and  a  posteriori  error  estimates  for  1-D,  2-D  and  3-D  BVP 
and  IVP  that  are  dependent  on  independent  computational  parameters  h,  p  and  k. 

The  work  is  basic  and  fundamental  in  that  it:  (i)  requires  a  closer  examination  of  the 
global  behaviors  of  approximations  to  assess  and  determine  the  intrinsically  important  ap¬ 
proximation  features  dependent  on  k  that  are  ignored  in  the  currently  used  methodologies, 
(ii)  requires  development  of  fundamental  approaches  through  with  approximation  aspects 
dependent  on  k  can  be  incorporated  in  the  a  priori  and  a  posteriori  error  estimates.  Our 
preliminary  investigations  show  that  inter-subdomain  behaviors  of  the  approximations 
(normal  derivatives  of  dependent  variables)  ignored  in  the  currently  derived  estimates 
are  dependent  on  k.  These  behaviors  must  be  incorporated  in  the  error  estimates  of  the 
quantities  of  interest  so  that  the  resulting  estimates  will  hold  globally  in  the  pointwise 
sense. 

A  priori  and  a  posteriori  estimates  for  BVP  and  IVP  must  be  considered  in  h,  p,  k 
framework  in  which  approximations  are  in  higher  order  scalar  product  spaces  Hk,p  (k 
being  the  order  of  the  space).  For  initial  value  problems,  such  estimates  must  consider 
k  =  (ki,  k2),p  =  (pi > P2)  in  which  pi  >  2 k\  —  l,p2  >  2A;2  —  1  in  which  k\  and  k2  are  the 
orders  of  the  approximation  spaces  in  space  and  time  and  pi ,  p2  are  the  degrees  of  local 
approximations  in  space  and  time.  The  importance  of  these  estimates  in  adaptive  pro¬ 
cess,  estimation  of  theoretical  convergence  rates  achievable,  and  rigor  of  the  mathematical 
framework  for  finite  element  processes  has  been  well  recognized.  The  concepts  discussed 
here  will  bring  the  mathematical  foundation  for  finite  element  computations  in  h,  p,  k 
framework  to  a  more  complete  stage  in  which  adaptive  processes  in  h,  p,  k  framework  can 
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be  initiated  and  designed  for  optimality  of  computational  efficiency  and  appropriate  error 
control. 


3.2  A  review  of  current  methodology  and  development  of  con¬ 
cepts  in  hpk  framework 

In  this  section  we  present  some  preliminaries  and  a  short  review  of  currently  used  method¬ 
ology  of  deriving  the  estimates  and  demonstrate  the  need  for  new  research  in  this  area. 
The  finite  element  approximations  in  Sobolev  space  (generally  of  order  1)  and  the  asso¬ 
ciated  infrastructure  leading  to  the  theory  of  generalized  solutions  or  theory  of  distribu¬ 
tions  forms  the  mathematical  backbone  of  a  priori  and  posteriori  error  estimations.  To 
demonstrate  the  features  of  the  currently  used  mathematical  framework,  we  consider  the 
following  1-D  simple  boundary  value  problem, 


~i£  =  f  in  n  =  |0-i| 


with 


V?|0  =  yo  and  —  \L  =  q 


dtp 

dx 


Galerkin  method  with  weak  form  yields  the  following  integral  form, 


(1.1) 

(1.2) 


B{v,ip)n  =  l{v)n  (1.3) 

in  which 

B{v’v)n  =  Ltxfxia  (L4) 

l(v)n  ~  (f,v)n  +  qv(L)  (1.5) 

and  v  —  Sip;  () 

Let  fr  =  (Uefje)  be  a  discretization  of  O  in  which  f le  is  a  typical  subdomain  ’e’  of  uT 
and  let  ifh  =  U eipeh  be  an  approximation  of  ip  in  BlF  where  ip\  approximates  p>h  over  an 
element  ’e’.  The  ipn  would  be  a  solution  of 


B{y,iph)  =  l(v)  ;  v  =  6iph  (1.6) 

For  the  discretization  QT ,  we  obtain 

a  =  E'eW  ;  «  =  Svl  (1.7) 

e  e 

In  (1.7),  Be(v,  iph)  is  obviously  from  (1.4).  However,  le[y)  requires  a  closer  examination. 
Using  (1.5),  we  obtain 

£  m  =  £(/,  v)a.  +  £(»M)Pf  +  v(4 )/%))  (14) 

e  e  e 


where  x\  and  x\  are  nodal  coordinates  of  the  two  end  nodes  of  the  element. 
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Remarks 

(1)  If  </?/(  and  hence  ipeh  are  of  class  C°,  i.e.  in  HliP(Qe)  scalar  product  space,  then  ^ 
is  discontinuous  at  the  inter-element  boundaries  and  (1.7)  from  (1.3)  is  not  possible 
in  the  strict  sense  of  calculus  of  continuous  and  differentiable  functions. 

(2)  However,  based  on  theory  of  distributions,  (1.7)  from  (1.3)  is  valid  in  jy1,p(F!e)  space 
in  which  (ph.  is  of  class  C°.  In  such  approach  it  is  well  known  that  inter-element 
behavior  is  ignored. 

(3)  Furthermore,  when  computing  for  nodal  degrees  of  freedom  in  iph  (or  peh  i.e.  con¬ 
stants  in  the  linear  combination)  using  (1.7),  one  must  impose  the  condition  that 
sum  of  secondary  variables  at  the  inter-element  boundaries  is  zero  (assuming  no  ex¬ 
ternal  disturbance  is  applied  at  such  nodes).  This  condition  is  not  supported  by  the 
local  approximations  <peh  over  Qe  when  peh  are  of  class  C°,  due  to  the  fact  that  ^ 
is  discontinuous  at  the  interelement  boundaries.  Thus  when  < peh  are  of  class  C°,  we 
have  a  computational  process  in  which  there  are  inconsistencies  at  all  inter-element 
boundaries.  If  one  involves  h,  p,  fop-adaptive  processes  with  p>\  of  class  C°,  then  a 
vast  amount  of  resources  are  spent  in  overcoming  the  inconsistency  described  here. 

(4)  When  <p%  are  of  class  C1,  i.e.  in  7/2,p(Qe)  scalar  product  space:  (a)  (1.7)  from  (1.3) 
conforms  to  the  calculus  of  continuous  and  differentiable  functions  due  to  the  fact 
that  (fih  and  ^  are  continuous  everywhere  in  including  inter-element  boundaries, 
and  hence  the  integrand  in  (1.3)  is  continuous  in  ClT  (b)  the  conditions  on  the  sum 
of  the  secondary  variables  is  enforced  correctly  at  the  inter-element  nodes  due  to  the 
fact  that  local  approximations  peh  in  //2,p(fie)  space  support  this  condition.  Hence, 
there  is  no  inconsistency  in  the  computational  process  when  are  of  class  C2  in 
H2’p{ne)  space.  The  issues  discussed  here  become  far  more  serious  and  damaging 
in  two  and  three  dimensions 

This  simple  model  BVP,  its  finite  element  formulation  and  the  remarks  presented 
here  demonstrate  how  one  ignores  the  inter-element  behaviors  and  creates  inconsis¬ 
tency  in  the  resulting  computational  process  in  one  approach  (current  approach)  and 
how  simply  these  are  avoided  all  together  in  the  hpk  framework.  Since,  the  a  priori 
and  posteriori  error  estimation  techniques  are  based  on  weak  forms  such  as  (1.6)  or 
(1.7)  in  which  are  of  class  C°  in  Hl,p{ Oe)  spaces,  it  is  rather  straightforward  to 
observe  that  in  these  methodologies  inter-element  behaviors  are  ignored  (just  like  it 
is  in  (1.6)  or  (1.7)  where  <peh  is  of  class  C°).  A  priori  and  a  posteriori  error  estimates 
used  currently  are  derived  in  terms  of  Lebesgue  measures  that  do  not  incorporate 
the  inter-element  behaviours.  Such  estimates  thus  naturally  are  only  functions  of 
h,p  and  smoothness  of  the  theoretical  solution  in  each  element.  In  view  of  the  fact 
that  all  finite  clement  computations  are  actually  dependent  on  /?.,  p  and  k,  the  need 
for  the  new  research  in  h,  p,  k  framework  is  quite  clear. 

(5)  Our  preliminary  work  shows  that  inclusion  of  inter-subdomain  behaviors  in  the  a  pri¬ 
ori  and  posteriori  error  estimates  is  important  and  crucial.  This  can  be  demonstrated 
by  considering  a  simple  1-D  case  (like  described  above)  in  which  the  discretization 
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is  uniform  with  characteristic  length  h  and  the  degree  of  local  approximation  is 
p.  Let  the  local  approximation  be  of  class  C°.  We  wish  to  establish  a  bound  on 
||  ^  || l2.  In  this  case  ^  exhibits  discontinuity  (jump)  at  the  inter-element 

boundaries.  One  finds  that  the  rate  of  convergence  of  the  error  in  the  jumps  of  ^ 
is  lower  than  the  rate  of  convergence  of  the  errors  in  ^  over  the  interiors  of  the 
subdomains.  The  same  holds  true  for  other  measures  as  well.  Thus  the  measures 
based  strictly  on  the  open  domains  of  the  elements  are  in  fact  over  estimates  of  the 
local  features  of  the  solution  that  are  often  of  great  interest  to  the  analyst  (such  as 
stresses  or  fluxes  across  the  interfaces) 

(6)  A  clear  demonstration  of  how  the  inter-element  behaviors  of  the  solution  derivatives 

is  ignored  in  the  presently  used  estimates  is  presented  in  this  section.  For  simplicity, 
consider  a  one-dimensional  problem  similar  to  one  described  earlier.  Let  h  and  p  be 
the  characteristic  length  and  degree  of  local  approximation.  Let  <p/(  be  of  class  C°. 
Then,  for  each  fle  the  estimates  °f  II  ^  \\l2  only  depend  upon  h  and  p.  Hence 

estimates  of  ||  ^  ^  \\l2  over  QT  if  derived  solely  based  on  ||  ^  —  ^7  ||i2  over  fie, 

will  naturally  depend  on  h  and  p  only  as  well  as  the  smoothness  of  the  theoretical 
solution  ip.  In  this  process  inter-element  behaviors  of  (^  —  ^jr)  has  been  completely 
ignored.  Currently  derived  estimates  are  based  on  this  methodology. 

(7)  A  global  estimate  of  ||  ^  -  ^  ||l2  that  holds  true  for  entire  DT  requires  that  we 

incorporate  measures  of  (^  —  ^-)  at  the  inter-element  boundaries  in  ||  ^  |U2 

i.e.  the  estimates  must  be  derived  over  fP.  Inter-element  behaviors  of  the  solution 
derivatives  normal  to  the  inter-element  boundaries  is  a  function  of  h,p  as  well  as  k. 

If  the  local  approximations  are  of  class  C°  (k= 1),  then  ^  (i.e.  derivative  of  order 
k—1)  is  discontinuous  at  the  inter-element  boundaries.  When  ipeh  are  of  class  C1, 
then  (i.e.  derivatives  of  order  k= 2)  exhibit  discontinuity  at  the  inter-element 
boundaries.  That  is  a  measure  of  (^£  -  ^r)  at  the  inter-element  boundaries  is  a 
function  of  k  (as  well  as  h  and  p),  the  order  of  the  approximation  space. 

Thus,  a  priori  and  posteriori  error  estimates  must  be  in  h.,p,  k  framework  in  which 
the  estimates  will  incorporate  inter-element  behavior  and  thereby  including  k,  the 
order  of  the  approximation  space  in  the  estimates  in  addition  to  h  and  p. 

4  Rediscretizations,  moving  meshes  and  solution  map¬ 
ping  strategies  and  associated  computational  in¬ 
frastructure  for  BVPs  and  IVPs  in  h.pk  framework 

Appendix  A  contains  a  research  report  issued  under  ’Computational  Mechanics  Labora¬ 
tory’  of  the  Department  of  Mechanical  Engineering  of  the  University  of  Kansas,  Lawrence, 
KS.  This  report  describes  the  research  work  done  in  this  area  during  the  third  year  of  the 
grants. 
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5  Future  work 


The  fluid-solid  interaction  and  error  estimation  work  presented  here  is  preliminary.  At 
present,  many  Ph.D  students  are  engaged  in  utilizing  the  concepts  presented  here  to 
develop  a  complete  infrastructure  in  these  two  areas  in  which  ID,  2D  and  3D  BVPs 
and  IVPs  will  be  treated  rigorously  without  bias  to  their  origin  or  fields  of  application. 
hpk  adaptive  processes  are  also  currently  an  area  of  focus.  This  research  will  utilize  the 
developments  in  a  priori  and  a  posteriori  error  estimation  for  developing  methodologies 
and  algorithms  for  h,  p,  k  adaptive  processes. 

The  research  work  presented  here  on  rediscretizations,  moving  meshes  and  solution 
mapping  is  of  critical  significance  in  solid  mechanics  areas  of  large  deformation,  large 
strain  utilizing  Lagrangian  mathematical  models.  In  such  cases  discretizations  become 
excessively  distorted  during  evolution  (IVPs)  or  load  incrementing  (BVPs)  and  hence  it 
becomes  necessary  to  rediscretize  and  obtain  a  map  of  the  existing  solution  onto  the  re¬ 
discretization.  Moving  mesh  approach  permits  efficient  computations  in  which  moving 
micro  fronts  can  be  resolved  in  a  macro  domain.  The  future  work  in  this  area  consists 
of  some  algorithm  developments  but  largely  of  applications  in  various  areas  of  continuum 
mechanics  that  may  necessitate  specific  developments  related  to  those  areas.  As  an  exam¬ 
ple,  mapping  strategies  for  fracture,  plasticity  and  damage  mechanics  in  general  in  which 
the  quantities  of  interest  are  generally  accumulated  on  incremental  basis  are  generally 
history  dependent  and  hence  may  require  specific  developments. 
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ABSTRACT 


Rediscrctization,  moving  mesh  and  solution  mapping  approaches  are  presented  for  boundary  value 
problems  (BVPs)  and  initial  value  problems  (IVPs).  When  BVPs  are  solved  using  mathematical 
models  in  Lagrangian  frame  of  reference,  the  discretizations  often  become  excessively  distorted 
during  load  incrementing.  Rediscrctization  and  a  total  map  of  the  solution  from  the  original 
discretization  is  beneficial  in  such  cases.  In  the  case  of  IVPs,  using  space-time  strips  with  time 
marching,  the  moving  mesh  approach  permits  accurate  evolutions  of  high  solution  gradient  moving 
fronts  using  computationally  practical  spatial  discretizations.  Both  the  rediscrctization  and  moving 
mesh  procedure  require  mapping  of  the  geometry  of  the  new  discretization  (whole  or  part  of  the 
boundary)  onto  the  original  discretization  and  then  mapping  of  the  solution  using  this  geometry 
map. 

Rediscrctization  and  moving  mesh  approaches  for  ID  and  2D  cases  are  considered  for  local 
approximation  of  Lagrange  type,  p- version  hierarchical  as  well  as  C\ClJ  ap¬ 

proximations  of  higher  order  global  differentiability.  Merits  of  higher  order  global  differentiability 
approximations  in  h,p,k  mathematical  and  computational  framework  are  demonstrated  for  re- 
discretizations  as  well  as  moving  meshes.  Inter-element  flux  problems  present  in  currently  used 
procedures  arc  examined  for  C°  local  approximations.  Remedies  arc  presented  and  their  absence  is 
clearly  shown  when  the  local  approximations  arc  of  higher  classes  than  C°.  Introduction,  literature 
review  and  scope  of  work  is  presented  in  Chapter  1.  Chapter  2  contains  developments  of  method¬ 
ologies,  approaches  and  procedures.  Numerical  studies  are  presented  in  Chapter  3.  Summary  and 
conclusions  are  given  in  Chapter  4. 
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Chapter  1 


Introduction  and  Proposed 


Research 


Moving  mesh  and  rcdiscretization  strategics  are  important  and  an  integral  part  of  practical  and 
accurate  computational  methodologies.  The  need  for  such  approaches  arise  in  boundary  value 
problems  (BVP)  as  well  as  initial  value  problems  (IVP).  These  strategies  must  be  considered  and 
evaluated  for  mathematical  models  constructed  in  Lagrangian  frame  of  reference  as  well  as  Eulerian 
frame  of  reference. 

In  Lagrangian  frame  of  reference  the  reference  frame  and  the  observer  both  are  stationary 
and  all  quantities  of  interest  are  measured  in  the  fixed  frame  of  reference.  In  this  approach,  the 
material  particles  and  the  discretization  are  the  same  and  hence  as  the  material  particles  undergo 
large  motion  and  large  strain  the  discretization  become  progressively  distorted  due  to  which  the 
computations  deteriorate  and  may  eventually  cease.  In  order  for  the  computations  to  proceed  and 
to  be  reliable  and  accurate,  rediscretization  is  often  necessary  when  the  mesh  distortion  becomes 
excessive  and  the  current  state  of  the  solution  must  be  mapped  accurately  onto  the  rediscretization. 
The  mathematical  models  in  Lagrangian  frame  of  reference  are  most  commonly  used  for  solid 
continua. 

When  the  mathematical  models  are  constructed  in  Eulerian  frame  of  reference  and  simulated 
numerically,  the  discretization  remains  fixed  and  the  material  particles  flow  through  it.  This  ap¬ 
proach  is  ideally  suited  for  arbitrary  large  motion  and  large  strain  rates  and  hence  is  the  preferred 
strategy  of  constructing  mathematical  models  for  fluid  flow,  gas  dynamics,  polymer  flow,  etc.  The 
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most  significant  advantage  of  this  approach  is  that,  since  the  mesh  remains  fixed,  there  is  no  mesh 
distortion.  Hence,  in  the  case  of  BVPs,  there  is  very  little  need  for  rediscretizations  if  any.  How¬ 
ever,  in  many  instances  IVPs  contain  high  localized  gradient  moving  fronts  such  as  in  the  case  of 
Riemann  shock  tube,  combustion  process  with  flame  propagation,  etc.  In  such  cases,  severe  mesh 
grading  may  be  required  in  the  zones  containing  the  fronts  but  a  relatively  coarser  discretization 
can  be  employed  elsewhere.  Since  the  fronts  propagate  as  time  evolution  proceeds,  the  high  fidelity 
of  the  solution  in  the  local  zone  containing  fronts  necessitates  the  same  level  of  discretization  ahead 
of  the  fronts.  This  often  results  in  discretizations  that  are  computationally  impractical.  For  such 
applications  a  space-time  time  marching  approach  is  highly  meritorious  [1],  In  this  approach  one 
constructs  a  space  time  strip  or  a  slab  for  an  increment  of  time  and  computes  a  converged  space¬ 
time  solution  for  it  through  h,p,k  adaptivity.  The  solution  state  at  the  open  boundary  (t  =  to  +  At), 
to  being  initial  time)  serves  as  the  initial  condition  for  the  second  space-time  strip  or  slab  for  which 
a  converged  solution  is  obtained.  This  process  is  continued  until  the  desired  time  is  reached.  In 
this  strategy  it  is  obvious  that  the  discretization  in  space  for  the  second  space-time  strip  needs  to 
be  of  the  same  level  of  refinement  as  that  for  the  time  t0  containing  the  front  if  high  fidelity  of  the 
fronts  is  to  be  maintained.  The  same  holds  true  for  subsequent  stages  of  the  evolution. 

An  alternative  to  the  uniformly  fine  discretization  in  the  zones  containing  fronts  and  ahead 
of  the  front  is  to  move  the  discretization  at  time  to  with  the  same  speed  as  the  speed  at  which 
the  fronts  are  moving.  This  automatically  provides  the  needed  refinement  for  the  front  to  evolve 
accurately  during  the  subsequent  time  steps.  An  advantage  to  this  approach  is  that  the  total 
degrees  of  freedom  (dofs)  at  each  space-time  strip  or  slab  remain  approximately  the  same. 

From  the  standpoint  of  mathematical  treatment  the  rcdiscretization  and  moving  mesh  ap¬ 
proaches  are  almost  identical.  In  both  cases  the  following  features  are  common. 

(1)  There  is  an  existing  discretization  with  numerical  results,  i.e.  numerical  values  of  all  dofs  and 
the  complete  solution  state  arc  known. 

(2)  There  is  a  new  discretization  which  may  be  arbitrary  (for  the  sake  of  generality)  compared 
to  the  existing  discretization  for  which  the  numerical  solution  is  not  known. 
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(3)  The  geometry  of  the  new  discretization  needs  to  be  mapped  on  to  the  existing  discretization 

(4)  The  complete  state  of  the  numerical  solution  needs  to  be  mapped  from  the  existing  discretiza¬ 
tion  onto  the  new  discretization. 

(5)  In  the  case  of  BVPs  a  total  map  is  required  for  the  entire  new  discretization. 

(6)  In  the  case  of  IVPs  with  a  space-time  strip  or  slab  and  space-time  time  marching,  only  the 
solution  at  the  open  boundary  of  the  previous  space-time  strip  needs  to  be  mapped  as  the 
initial  condition  for  the  new  space-time  strip. 

(7)  Development  of  measures  to  ensure  that  the  solution  mapping  from  the  existing  discretization 
onto  the  new  discretization  accurate. 

In  principle,  the  rcdiscrctization  and  moving  mesh  methodologies  appear  rather  simple  and  straight¬ 
forward.  However,  their  success  depends  upon  the  accuracy  of  the  solution  mapping. 


1.1  Basic  Concepts,  Methodologies  and  Literature  Review 

Functional  analysis  in  Sobolev  spaces  and  approximation  theory  basically  constitute  the  mathe¬ 
matical  foundation  of  the  currently  used  finite  element  method.  In  this  approach,  local  approxima¬ 
tions  over  subdomains  (elements)  of  class  C°  and  Galerkin  method  with  weak  form  have  primarily 
dominated  the  research  over  the  past  forty  years.  The  problems  in  the  computational  processes 
associated  with  Galerkin  method  with  weak  form  when  the  differential  operators  (for  both  BVP 
and  IVP)  arc  non-self-adjoint  and  non-linear  have  been  reported  by  Surana  et.  al.  [2-4]  and  arc 
not  a  subject  of  study  in  the  present  work.  However,  the  role  of  C°  local  approximations  in  rcdis¬ 
crctization  and  moving  mesh  strategics  is  of  critical  importance  and  its  consequences  need  to  be 
examined  in  detail. 

First,  we  note  that  when  the  local  approximations  are  of  class  C°,  the  global  approximation 
for  the  whole  discretization  is  also  of  the  class  C°.  Such  local  approximations  exhibit  discontinuity 
in  the  first  derivatives  of  the  dependent  variables  normal  to  the  inter-element  boundaries  which 
cannot  be  quantified  in  the  sense  of  Z/2-norms  due  to  the  fact  that  such  discontinuities  occur 
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on  sets  of  measure  zero.  Instead,  their  influence  can  only  be  measured  in  terms  of  how  they 
affect  the  Z/2-norrns  over  the  interiors  of  the  subdomain.  This  is  referred  to  as  inter-element  flux 
discontinuity  problem.  When  solutions  with  such  inter-element  flux  discontinuities  are  mapped 
onto  a  new  discretization  in  which  the  inter-element  boundaries  do  not  coincide  with  those  in 
the  existing  discretization,  new  inter-element  flux  discontinuities  are  created  at  different  locations. 
Furthermore,  the  jump  discontinuity  magnitude  in  the  map  may  be  quite  different  than  those  in 
the  original  discretization.  This  problem  is  well  recognized  [5]  and  there  are  many  published  work 
that  present  various  procedure  to  minimize  such  problems.  Nonetheless,  it  is  well  recognized  and 
agreed  upon  that  such  flux  problems  always  exist  in  almost  all  rediscretization  and  moving  mesh 
procedures.  Some  pertinent  published  work  is  discussed  in  the  following. 

As  described  above,  a  need  for  rcdiscretization  or  moving  mesh  method  typically  arises  when 
solving  BVPs  with  a  Lagrangian  description  or  IVPs  with  an  Eulcrian  description.  A  purely 
Lagrangian  approach  is  limited  by  its  ability  to  maintain  reliable  computations  when  the  material 
and  hence,  the  mesh  undergo  large  distortion.  It  has  been  suggested  by  Hughes  et.  al  [5]  to  combine 
the  best  features  of  both  approaches  in  the  Arbitrary-Lagrangian-Eulerian  (ALE)  method.  The 
ALE  method  employs  continuous  rezoning  of  the  mesh  to  provide  both  flexibility  for  dealing  with 
large  distortions  in  the  material  and  resolution  for  capturing  details  in  the  solution.  This  proves 
especially  useful  when  applied  to  fluid-structure  interaction  problems.  The  difficulty  in  the  ALE 
method  is  the  difficulty  in  finding  algorithms  to  prescribe  the  appropriate  nodal  displacements 
during  rezoning.  The  authors  have  noted  the  need  for  studies  regarding  the  accuracy  for  such  an 
algorithm. 

The  Moving  Finite  Elements  (MFE)  of  Miller  [6,  7]  has  been  developed  to  maintain  resolution 
for  IVPs  with  sharp  solution  gradients  or  ”  near-shocks” .  In  such  generalized  methods,  nodes 
automatically  move  to  critical  areas  where  finer  resolution  is  needed  to  accurately  describe  the 
solution.  The  usefulness  of  MFE  has  been  demonstrated  in  the  solution  of  transient  Burgers 
equation  in  which  the  nodes  rapidly  accumulate  where  the  solution  is  steepest.  A  concern  of  MFE 
is  the  possibility  of  overturned  solutions  or  shocks,  as  discussed  for  convection  diffusion  phenomena 
by  Baines  [8].  Baines  has  reported  great  success  in  solving  difficult  problems  when  employing 
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appropriate  penalty  functions.  The  inter-element  flux  problems  inherent  in  these  procedures  due 
to  C°  local  approximations  remain  unresolved. 

1.2  Scope  of  Present  Research 

In  the  present  work,  rcdiscretization  and  moving  mesh  procedures  are  developed,  evaluated  and 
numerically  demonstrated  in  h,p,k  mathematical  and  computational  framework  [1-4],  In  this  frame¬ 
work,  the  order  k  of  the  approximation  space  is  an  independent  computational  parameter  in  addition 
to  h,  the  characteristic  discretization  length  and  p,  the  degree  of  local  approximation.  The  param¬ 
eter  k  defines  global  differentiability  of  the  approximation.  Thus,  the  local  approximations  are 
considered  in  Hk'p{ f2e)  scalar  product  spaces  containing  approximation  functions  of  degree  p  yield¬ 
ing  global  differentiability  of  order  ( k  -  1),  i.e.  global  approximations  of  class  Ck.  Methodologies 
and  procedure  are  developed  (Chapter  2)  for:  (1)  mapping  the  new  discretization  onto  the  existing 
discretization  (2)  mapping  of  the  solution  from  the  existing  discretization  onto  the  rediscrctization 
(3)  measures  that  ensure  the  accuracy  of  the  mapped  solution  for  ID  and  2D  BVPs  as  well  as 
IVPs  in  h,p,k  framework.  Numerical  studies  for  ID  and  2D  BVPs  and  ID  IVPs  are  presented  in 
Chapter  3.  A  summary  of  the  work,  conclusions  and  some  directions  for  future  work  are  presented 
in  Chapter  4. 
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Chapter  2 


Rediscretization  and  Moving  Mesh 
Methodology:  Mathematical  and 
Computational  Developments 


This  chapter  contains  development  of  mathematical  and  computational  methodologies  for  redis- 
crctization  and  moving  meshes  as  well  as  developments  of  the  measures  to  assess  the  validity  and 
accuracy  of  the  proposed  approach.  First,  we  consider  rcdiscrctizations  associated  with  BVPs.  For 
the  sake  of  simplicity,  we  consider  a  BVP  in  one  dependent  and  independent  variable  (an  ordinary 
differential  equation).  One-dimensional  axial  deformation  of  an  elastic  rod,  convection-diffusion 
equation  and  Burgers  equation  are  a  few  examples  of  such  BVPs. 

2.1  Approximations  of  Class  C° ,  Lagrange  Functions,  ID  Case 

Let 

Acf>-f  =  0  in  P=(0  ,L)  (2.1) 

be  the  BVP  (with  some  boundary  conditions).  Let  the  differential  operator  be  of  order  2m  and  also 
let  the  theoretical  solution  of  (2.1)  be  of  class  Cfc(fi),  k  >2.  A  finite  element  formulation  of  (2.1) 
may  be  constructed  using  any  appropriate  integral  strategy  based  on  the  nature  of  the  differential 
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Figure  2.1:  Discretization  ml  (a):  current  mesh  and  rediscrctization  m2  (b):  new  mesh 
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operator  A  [2]- [4] .  Let  be  a  four-element  discretization:  mi  (Figure  2.1)  for  obtaining  a 
numerical  solution  of  (2.1),  in  which  is  a  subdomain  (or  an  element)  such  that 

^  =  (2’2) 
e 

Let  be  the  global  solution  (approximation)  for  D©  Then: 

=  (2.3) 

e 

in  which  ml4 ’hi&mi)  the  local  approximation  for  a  subdomain  S2®nl.  For  the  sake  of  simplicity 

we  assume  the  p-levcl  to  be  two  and  the  approximation  functions  to  be  of  standard  Lagrange  type 
implying  the  function  values  at  each  node  arc  the  nodal  degrees  of  freedom  and  ml4>h  and 
arc  of  class  Co-  Let  each  element  of  the  discretization  be  mapped  in  the  natural  coordinate  space 
f;  with  local  node  numbers  1,2  and  3  and  let  rnlxei\  i  =  1,2,3  to  be  the  global  coordinates  of  the 
nodes  of  the  clement  (Figure  2.2).  Let 

3 

mV(^)  =  £^i(£,,?r1X?  (2-4) 

i— 1 
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(a)  element  ©  of  fiTra,  in  (b)  map  of  element  ©  of 

physical  coordinate  space  £iTml  in  natural  coordinate 

space 

Figure  2.2:  Element  e  and  its  map  in  natural  coordinate  space 

define  the  mapping  of  element  e  of  0^  with  domain  in  the  natural  coordinate  space. 

Figure  2.1  (b)  shows  a  four  element  rcdiscretization  (m2)  of  the  domain  D  in  which  the  interior 
nodes  (including  intcr-elemcnt  boundaries)  do  not  coincide  with  the  original  discretization  mi  (for 
the  sake  of  generality). 

Let  be  the  map  of  ”,V/t(^mi)  over  the  discretization  m2  such  that 

m2^(^2)  =  Um2^^m2)  (2.5) 

e 

in  which  m20/,(n^2)  is  the  approximation  of  d  over  ,  an  element  e  of  the  discretization  m2. 
Determination  of  m20/l(D^2)  ovcr  &m2  from  ml  over  is  a  two  step  process.  First  we 

determine  the  locations  of  the  nodes  of  mesh  m2  onto  mesh  ml.  We  refer  to  this  step  as  mapping 
of  rcdiscretization  m2  onto  discretization  ml.  In  the  second  step  we  determine  m2^/((D^2)  from 
rnl using  the  mapped  geometry.  Details  arc  described  in  the  following. 

2.1.1  Mapping  of  Rediscretization  (m2)  onto  Original  Discretization  (ml) 

First,  wc  note  that  the  Cartesian  coordinates  of  nodes  of  rcdiscretization  m2  arc  known.  Secondly, 
our  ultimate  goal  is  the  determination  of  from  ml rf’hi&mi) ■  This  essentially  requires  the 

determination  of  the  numerical  values  of  the  nodal  degrees  of  freedom  for  each  node  of  m2.  Once 
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the  nodal  dofs  for  m2  are  known,  m20|(fi^,2)  is  defined  through  local  approximations  and  hence 
m2<j) hi&mi)  by  (2-5)-  We  proceed  as  follows. 

(a)  Consider  an  clement  e  of  with  domain  D^2  and  nodal  coordinates  m2x®. 

(b)  Determine  the  elements  in  discretization  ml  containing  these  nodes.  Let  m2x-\i  =  1,2,3  lie 
in  elements  p,  q  and  r  of  respectively. 

(c)  Determine  the  £  coordinates  in  elements  p,  q,r  for  the  locations  m2x?,i  —  1,2,3.  Using  the 
geometry  mapping  for  elements  p,q,r  of  we  have  the  following  in  which  m2xf,  mlxf> 
mlXi>  mlXi  arc  known  ^-coordinates. 
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m2Xi  =  £  JVi(£)mlx? 

i=  1 

(2.6) 

m2X2  =  £M(£)mlx* 

i-i 

(2.7) 

3 

m2xl  =  T,^mlxl 

(2.8) 

i= 1 


Each  of  the  (2.6)-(2.8)  are  quadratic  expressions  in  £  and  hence  we  can  determine  a  pair  of  £ 
values  from  each  of  (2.6)-(2.8).  Let  (ml£p)j,f  =  1,2;  (ml£9)iti  =  1,2;  and  (ml£r)i,2  =  1,2  be 
the  pair  of  values  of  £  determined  from  (2.6)-(2.8)  respectively.  Discarding  those  outside  the 
range  [—1, 1]  ,  wc  obtain  the  £  coordinates  (ml£j)e,  (ml£2)e  and  (ml£3)e  for  the  ^-locations. 

(d)  Using  the  procedure  described  in  (c),  wc  can  determine  the  elements  of  &Jnl  containing  each 
node  of  the  elements  of  D^2  and  their  corresponding  £  locations  in  the  elements  of 

2.1.2  Mapping  of  onto  ClJn2,  i.e.  Determination  of  m2(/)/t(f2^2) 

For  an  element  e  of  U^2  with  the  map  of  its  nodes  in  elements  p,  q  and  r  of  Dmi  and  their  £- 
coordinates  determined  in  (c)  of  Section  2.2.1,  we  proceed  as  follows.  For  elements  p,  q  and  r  of 
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the  local  approximations  are  given  by 


3 

ml4(^mi)  =  53^(0  ml^';  i  =  V^r  (1 2 3 4-9) 

i=l 

where  ml  </>'•;  .7  —  1, 2, 3  are  nodal  dofs  for  element  i  of  discretization  ml  (known  function  values  in 
this  case).  We  note  that 

ml4A)l*=r^r  =  m<24>\  (2.io) 

mlK^qmi)k=(^r  =  m24>e2  (2.ii) 

my,(jCi  )|c=(«i6).  =  m2^  (2.12) 

Using  (2.9)-(2.12)  numerical  values  of  the  nodal  dofs  m2(j) f,  m2(j) m2<\>\  of  element  e  of  discretization 
are  obtained.  This  procedure  is  repeated  for  each  element  and  thus  eventually  all  numerical 
values  of  the  dofs  for  each  node  of  the  discretization  ClJn2  are  determined.  Now  using  local  approx¬ 
imations  for  each  element  the  solution  m2<j)eh  is  established  and  we  have 

3 

m24(^m2)  =  E^^)m2^  (2-13) 

i=l 

and 

m2<t»A)  =  U  m24fc)  (2.14) 

e 

Thus  we  have  the  map  ont°  from  mVft(^m l)  ont°  &mi- 

Remarks 

(1)  Many  issues  need  to  be  addressed  regarding  the  validity  and  accuracy  of  the  solution  mapping 
procedure  described  above. 

(2)  Extension  of  this  procedure  to  the  C°  p- version  hierarchical  local  approximations  is  needed. 

(3)  Extension  to  higher  order  continuity  local  approximations  is  required. 

(4)  Measures  to  assess  the  accuracy  of  the  mapping  procedure  need  to  be  established. 


10 


(5)  Extensions  to  2D  and  3D. 


2.1.3  Measures  of  the  Accuracy  of  the  Mapping 

In  the  development  presented  here,  there  arc  some  specific,  points  to  be  noted  that  may  be  helpful 
in  further  understanding  the  behaviors  of  mV/i  and  m2<f>h  and  development  of  the  measures  of 
accuracy  of  m2<j>h- 

(1)  and  m2^(D^2)  arc  of  class  C 2  due  to  p-lcvel  of  2,  however  ml  <f>eh(&mi)  and 

m20^t(D^l2)  are  °f  class  C°.  Thus,  and  ml eph {&mi)  are  °f  class  C2  over  = 

and  °m2  =  int{ Ue  fi^2)  but  of  class  C°  at  the  inter-element  boundaries. 

(2)  Since  and  rn24)PjL  are  of  class  C°  globally,  inter-clement  discontinuities  (jumps)  exist  in 
their  first  derivatives  normal  to  the  inter-element  boundaries  in  both  discretizations.  Loca¬ 
tions  of  the  jump  discontinuities  are  different  in  ml  and  m2.  The  magnitudes  of  the  jump 
discontinuities  may  also  differ  between  ml  and  m2. 

(3)  Based  on  (1)  and  (2),  d  ^ ,  ?h ^  and  d  -v-t'*-  ;  i  =  0, 1,  in  general  may  be  the  first 

ax  i2 

good  indicators  regarding  the  reliability  of  mapping.  Their  agreement  between  ml  and  m2 
discretizations  indicate  a  good  mapping.  Lack  of  agreement,  on  the  other  hand,  points  out 
the  existence  of  flux  discontinuity  problems. 

(4)  When  to  map:  (3)  raises  a  very  significant  question,  i.e.  when  to  map.  We  elaborate  on  this 
point  in  the  following  and  provide  some  clear  guidelines. 


(a)  Since  in  the  above  discussion  the  approximations  are  of  class  C°,  the  inter-element  jumps 


in  the  first  derivatives  exist.  However,  if  the  theoretical  solution  is  sufficiently  smooth  we 
expect  weak  convergence  of  the  class  C°  to  C 1  upon  h,p  refinements.  In  such  converged 


solutions,  the  jump  discontinuities  are  no  longer  significant  hence  we  expect 


vr-Uh) 

dx% 


dxl 


;  i  =  1, 2,  to  be  in  relatively  good  agreement. 


(b)  However,  when  is  not  weakly  converged,  jump  discontinuities  at  the  inter-element 
boundaries  arc  significant.  In  such  case,  mapping  the  solution  to  rediscretizations  is  not 


a  preferred  strategy. 
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(c)  Based  on  (a)  and  (b),  it  is  rather  straight-forward  to  sec  that  when 


(2.15) 


the  map  of  ml4>h  to  m2<Pk  is  unique  and  so  are  their  corresponding  Z/2-norms.  Likewise, 
when 


d(mlcph) 

d<p 

dx 

dx 

(2.16) 


the  maps  of  d ^  to 


dr2M 

dx 


arc  unique  and  so  arc  their  corresponding  L2-norms. 


(d)  It  is  rather  clear  from  (b)  and  (c)  that  before  attempting  to  map  the  solution  ml<f>h  on 
fi^2,  it  is  preferable  to  achieve  weak  convergence  of  mi(j>h  to  <p.  This  eliminates  the 
mapping  errors  due  to  the  fact  that  ml  0/,  =  <j).  Weakly  converged  solutions  for  higher 
p- level  Lagrange  elements  (with  function  values  as  nodal  dofs)  shall  permit  good  maps 
onto  f^2  in  which  higher  order  global  differentiability  is  possible  to  preserve  in  the  weak 
sense. 


2.2  p-version  Hierarchical  Local  Approximations  of  Class  C°:  ID 
case 

Consider  discretization  ml  of  Figure  2.1  in  which  the  elements  arc  three  node  p- version  C°  hierar¬ 
chical  elements  with  p- level  of  p^.  In  this  case  is  of  class  p^  over  f^r[1  =  *nt{ Ue  ^ml)  but 

is  of  class  C°  at  the  inter-element  boundaries.  Mapping  of  S2^i2  onto  follows  in  the  same  pro¬ 
cedure  described  in  Section  2.1.1.  However,  mapping  of  ml0ft(^mi)  onto  ^m2  to  obtain  m20/l(H^2) 
is  somewhat  different  than  the  procedure  discussed  in  section  2.1.2.  We  note  that  the  dofs  at 
the  end  nodes  are  function  values  whereas  at  the  mid-side  node  of  each  element  the  dofs  arc 
i  =  2, . . .  ,p£.  Thus,  at  each  mid-side  node  of  elements  0.fm2  numerical  values  of  needs  to  be 
obtained  from  the  elements  of  discretization  This  cannot  be  done  accurately  due  to  the  fact 
that  is  °f  class  C'°globally. 

This  situation  is  not  as  hopeless  as  it  might  appear.  A  careful  examination  of  the  derivation 
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Figure  2.3:  C°  p-version  (a)  non-hierarchical  and  (b)  hierarchical  Lagrange  elements 

of  C°  p- version  hierarchical  basis  functions  leads  to  a  resolution.  Figure  2.3  shows  C°  p- version 
non-hierarchical  Lagrange  elements  and  the  corresponding  p- version  hierarchical  Lagrange  elements 
derived  from  them  for  various  p- levels.  Let  us  assume  that  mVh(^mi)  's  a  converged  C°  solution 
for  discretization  ml  at  p- levels  of  p^.  It  is  desired  to  map  this  solution  onto  rediscretization  m2. 
This  obviously  requires  nodal  values  of  (j>  at  the  non-hierarchical  nodes  and  i  —  2 , ,p^  at 
each  of  the  hierarchical  nodes.  Since  the  C°  p-version  hierarchical  approximation  functions  are 
derived  from  the  corresponding  non-hierarchical  Lagrange  element  with  (p^  +  1)  equally  spaced 
nodes  between  —  1  <  £  <  1,  we  can  proceed  as  follows. 

(a)  Consider  0%n2  of  0^l2  with  the  end  nodes  located  at  m2x\  and  m2X3- 

(b)  Determine  elements  p  and  q  of  containing  these  nodes  with  the  corresponding  values  of 
£,  i.e.  (m2^)e  and  (m2{$)e. 

(c)  Consider  (p^  +  1)  equally  spaced  points  between  m2xf  and  rn2X;i  and  determine  their  £  coor- 
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dinates  in  elements  p  and  q  of  fljnl 

(d)  Use  the  local  approximations  of  elements  p  and  q  of  and  the  £  coordinates  determined 
in  (c)  to  obtain  function  values  at  ( p ^  +  1)  equally  spaced  points. 

(e)  Convert  these  (p^  +  1)  function  values  to  4$;  i  =  2, . . .  ,p^,  at  the  m2X2  coordinate  location. 
This  completes  the  map  of  ml0(fi^1)  over  £lem2  °f  &m'i ■  This  procedure  is  repeated  for  each 
element  of  0^2  to  obtain 

Remarks 

(1)  When  ml (j)h(fljni)  of  class  C°  is  weakly  converged  to  good  accuracy,  ||d>  — ml  — >  0  and 

the  proximity  of  ml4>h  to  the  theoretical  solution  <j>  is  established. 

(2)  Based  on  (1),  computations  of  ~  for  the  hierarchical  nodes  of  0^2  using  function  values  at 
(p^  +  1)  equally  spaced  points  between  the  end  nodes  of  each  element  of  0^2  appear  to  be  a 
reasonable  proposition. 

(3)  In  this  procedure,  accuracy  of  the  function  is  the  best  in  the  map  and  it  deteriorates  for 
progressively  higher  order  derivative  hierarchical  dofs. 

(4)  This  procedure  is  expected  to  yield  better  accuracy  of  the  map  m2(f>h  as  compared  to  the 
procedure  in  which  ~  for  the  hierarchical  nodes  are  directly  computed  using  ml0/i(Umi)- 

(5)  Numerical  studies  related  to  this  procedure  are  not  presented  in  this  thesis  and  will  be  a 
subject  of  future  study. 

2.3  p-version  Approximations  of  Higher  Classes  (CJ;  j  =  1, . . .):  ID 
case 

Consider  discretizations  ml  and  m2  shown  in  Figure  2.1  in  which  the  local  approximations  for 
each  element  can  be  of  higher  class.  Consider  classes  C'7(P^1);  j  —  1,2,...  in  scalar 

product  spaces  Hk’p(fl^nl)  in  which  p  —  2k  -  l  and  j  =  k  —  1.  p-levels  p  =  2j  —  1  are  the  minimum 
p- levels  for  j  =  1, 2, . . ..  At  these  p-levcls  the  hierarchical  center  node  has  no  degrees  of  freedom. 
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Let  m2  be  a  rcdiscrctization  (Figure  2.1  (b))  in  which  the  locations  of  the  interior  nodes  of  the 
mesh  are  not  coincident  with  those  in  discretization  mi. 

Mapping  of  discretization  m2  onto  mi  follows  the  same  procedure  as  described  in  Section  2.1.1. 
For  mapping  of  the  known  solution  ml0/t(fimi)  onto  to  obtain  m24>h{ f^2)>  we  consider  the 
following  two  cases. 

(a)  Let  the  p- level  be  p  =  2j  —  1  where  j  is  the  order  of  continuity.  In  this  case,  only  the  end 
nodes  have  dofs.  We  note  that  since  mV/i(^mi)  °f  class  c3  in  f i ^  but  of  class  p  in  i.e. 
at  the  inter-clement  boundaries  the  solution  class  is  c3  but  in  the  interiors  of  the  element  the 
solution  class  is  cp.  When  is  mapped  onto  using  the  procedure  described  in  section 
2.1.2  wc  expect  c3  global  differentiability  features  to  be  preserved.  Consider  the  following  two 
cases. 

(i)  Discretization  ml  (i.e.  h  and  p)  is  such  that  ml<ph  is  not  a  converged  solution 

(ii)  Discretization  ml  is  refined  sufficiently  with  adequate  p-levels  to  obtain  a  converged 
solution  ml<t>h  before  mapping  it  onto  &Jn2. 

The  convergence  is  assumed  to  be  in  the  appropriate  norms  in  the  corresponding  approxima¬ 
tion  spaces.  Consequences  of  (i)  or  (ii)  are  discussed  and  illustrated  numerically  in  Chapter 
3. 

(b)  Consider  p-levels  p$  >  2j  +  1  where  j  is  the  order  of  continuity  of  ml</»/l(D^1).  In  the 
case  the  hierarchical  center  node  of  each  element  in  both  discretizations  has  dofs  that  are 

i  —  2j  —  1, . . .  ,p£.  We  consider  the  following  two  approaches  for  mapping  ml^/i(D^1) 
onto  n^2. 

(i)  Obtain  numerical  values  of  ^  for  the  hierarchical  nodes  of  Ujn2  directly  using  local 
approximations  for  the  elements  of  discretization  .  In  this  case  we  expect  the  same 
type  of  contamination  in  m20;,  as  described  for  C°  p- version  hierarchical  approximations 
even  when  the  solution  ml0/j  is  converged.  This  obviously  is  due  to  the  fact  that  global 
differentiability  of  ml(f>k  is  lower  than  what  is  required  for  obtaining  numerically  accurate 
values  of  hierarchical  dofs. 
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(ii)  Use  the  alternate  procedure  described  in  Section  2.1.2  which  employs  corresponding 
C°  Lagrange  configuration  to  obtain  In  this  approach  higher  accuracy  of  m2<j)h  is 
expected  due  to  the  fact  that  only  converged  function  values  are  utilized  in  the  compu¬ 
tations  of  hierarchical  dofs. 


2.4  Two-dimensional  Case:  Rediscretization 

Consider  discretizations  ml  and  m2  of  Fig.  2.4.  Let  mi  be  the  original  discretization  for  which 
the  solution  ml^(0^1)  is  known  and  m2  be  a  rediscretization  for  which  a  map  of 

is  desired.  Regardless  of  whether  the  local  approximation  of  type  C°  Lagrange,  C° 
p- version  hierarchical  or  that  of  higher  order  global  differentiability,  the  mapping  of  the  geometry 
of  rcdiscretization  m2  onto  ml  remain  the  same.  Hence,  we  consider  this  first. 


2.4.1  Mapping  of  Rediscretization  m2  onto  discretization  m\ 

Let  each  element  of  the  discretizations  ml  and  m2  be  mapped  onto  a  two-unit  square  natural 
coordinate  space  with  the  origin  of  the  coordinate  system  located  at  the  center  of  the  square. 
For  simplicity,  we  assume  the  shape  functions  in  the  mapping  to  be  quadratic.  For  a  typical  element 
e  of  njnl,  we  can  write 


mlxe 


(Cv)  = 


i=  1 
n 


mlxt 


mlyf 


(2.17) 


2=1 


To  determine  the  map  of  f onto  we  follow  the  procedure  described  in  Section  2.1.1  for  the 
ID  case.  The  basic  steps  arc  described  in  the  following. 

(a)  Consider  an  element  e  of  S7^2  with  domain  and  nodal  coordinates  (m2a;e,  m2xe)\  i  = 
1, . . . ,  n,  n  being  the  number  of  nodes. 


(b)  Determine  the  elements  in  the  discretization  containing  these  nodes.  Let  p,q,r,...  be 
such  elements. 
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(c)  It  is  sufficient  to  consider  a  typical  node  of  element  f2^2-  The  procedure  described  for  this 
node  can  be  repeated  for  the  remaining  node.  Consider  node  1  with  coordinates  (m2xf  ,m2  y\). 
Let  these  this  node  be  located  in  element  p  of  Determine  the  £,  77  coordinates  in  element 
p  for  the  location  (m2a;®,rn2  yf).  Using  the  geometry  mapping  (2.17)  for  element  p  of  we 
have 

71 


71—1 

n 

=  mi#f 


(2.18) 


71=  1 

in  which  (m2xe,m2  ye),  (mlx?™2  7/f);  i  =  l,...,n  are  known  coordinates.  Equations  (2.18) 
arc  quadratic  simultaneous  equations  in  £  and  77.  These  are  solved  numerically  using  Newton’s 
linear  method.  This  yields  two  roots,  i.e.  two  pairs  arc  f  and  i).  The  values  in  the  range 
[£,77]  =  [ —  1 , 1]  x  [—1,1]  are  retained  while  the  others  arc  discarded.  Let  ((ml£p)e,  (mlili)e) 
be  the  values  of  £  and  77  for  location  (m2xf,m2  yf)  in  element  p  of  x.  Using  this  procedure 
the  (£,77)  locations  for  all  nodes  of  f \em2  arc  determined  in  elements  p,q,r,...  of  Let 
((ml£i)e>  l)e)i  ((ml^I)e>  (ml??2)e)>  •  •  •  be  the  desired  ,  77  coordinates. 


(d)  This  procedure  described  in  (c)  is  repeated  for  each  element  of  0;^2,  thus  establishing  a  map 


of  onto  n£x. 


2.4.2  Mapping  of  mV/i(^mi)  onto  fi^2,  i.e.  Determination  of  Solution  Map 

It  is  sufficient  to  describe  the  solution  mapping  strategy  for  a  typical  element  e  of  tljn2  with  its 
domain  Slem2  and  the  £,  77  coordinates  of  its  map  in  elements  p,q,r,. . .  of  (determined  in  (c)  of 
Section  2.4.1).  For  a  typical  element  i  of  the  local  approximation  is  given  by 

nd 

=  (2.19) 

3= 1 

in  which  are  the  nd  dofs  for  element  i  of  Consider  the  following  types  of  local  approxi¬ 
mations  for  typical  elements  and  U^2  of  and  S7^2. 
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y 


Figure  2.4:  Discretization  ml  (a):  current  mesh  and  rediscretization  m2  (b):  new  mesh 


t  V  tn 


(a)  dement  ©  of  QTni]  in  (b)  map  of  element  ©  of 

physical  coordinate  space  QTml  in  natural  coordinate 

space 

Figure  2.5:  Element  e  and  its  map  in  natural  coordinate  space 
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(i)  <7°  local  approximations  of  Lagrange  type: 

Tensor  product  of  ID  (7°  Lagrange  functions  can  be  used  in  £,  r;  to  generate  the  approximation 
functions  for  such  elements.  For  such  elements,  each  node  of  every  element  has  only  function 
values  as  a  dof.  Hence,  using  local  approximations  of  the  type  (2.19)  for  the  elements  of  £ljnl 
and  the  £,  77  coordinates  of  the  nodes  of  the  elements  of  x,  the  function  values  can  be  easily 
determined.  This  procedure  is  exactly  parallel  to  that  described  for  ID  Lagrange  (7°  elements 
in  Section  2.1.2. 

(ii  )  C0-0  p-version  hierarchical  approximations: 

These  approximation  functions  can  be  derived  using  products  of  ID  p- version  hierarchical 
basis  functions  in  £  and  ?].  In  this  case  the  corner  nodes  are  non-hierarchical  nodes  and  have 
function  values  as  dofs.  Mid-side  and  center  nodes  on  the  other  hand  are  hierarchical  nodes 
and  contain  derivatives  of  the  function  with  respect  to  £  and  77  as  dofs  [9], [10].  We  can  proceed 
as  follows  and  possibly  develop  two  strategies. 

(a)  In  the  first  case,  dofs  at  the  non-hierarchical  nodes  are  determined  using  the  approach 

described  in  (i).  The  hierarchical  dofs  arc  also  determined  using  local  approximations 
for  the  corresponding  elements  of  The  problems  in  this  approach  are  exactly  the 
same  as  those  in  the  ID  case  in  Section  2.2.  It  appears  that  in  this  approach,  good 
maps  of  ml0/i(^ml)  onto  ^m2>  Le-  m2<Ph{^m 2),  are  not  possible.  That  is,  m2<ph{ *s 
always  contaminated  due  to  inaccurate  values  of  hierarchical  dofs,  which  essentially  is  a 
consequence  of  limited  global  differentiability  of  (<7°  in  this  case). 

(b)  Non-hierarchical  dofs,  i.e.  function  values,  are  determined  in  the  same  way  as  in  (a),  but 
hierarchical  dofs  arc  determined  using  their  corresponding  C°  Lagrange  configurations 
as  discussed  for  the  ID  case  in  Section  2.2.  This  approach  appears  meritorious  and  will 
be  investigated  further  in  future  work. 

(iii)  p-version  Approximations  of  Higher  Classes  (7ij;  i,j  =  1,2,...: 

In  this  category  wc  consider  two  different  types  of  local  approximations.  In  both  cases, 
element  nodal  configuration  is  that  of  standard  nine  nodes. 
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(a)  C‘ J  p- version  higher  order  continuity  approximations:  rectangular  elements 

If  the  element  shapes  in  the  physical  coordinate  space  arc  rectangular  and  if  the  coor¬ 
dinate  axes  x,y  are  parallel  to  £,  77,  then  the  approximation  functions  for  these  elements 
can  be  generated  using  tensor  products  of  ID  C\  type  p- version  approximations  in 
£  and  ?/  natural  coordinates  [10]. 

(b)  C,J  p- version  higher  order  continuity  approximations:  distorted  elements  in  the  physical 
coordinate  space 

When  the  element  shapes  in  the  physical  coordinate  space  are  distorted,  the  tensor  prod¬ 
uct  cannot  be  utilized  to  determine  approximation  functions.  In  this  case,  an  alternate 
procedure  [11]  is  necessary. 

Remarks 

(1)  The  dofs  at  the  corner  nodes  arc  function  values  as  well  as  derivatives  of  the  function  with 
respect  to  £  and  y  depending  upon  the  order  of  continuity  in  x  and  y. 

(2)  The  dofs  at  the  corner  nodes  for  elements  of  type  (a)  and  (b)  are  not  the  same,  the  major 
difference  being  that  all  dofs  at  the  corner  nodes  of  type  (a)  cannot  be  transformed  under 
pure  rotation  whereas  those  of  type  (b)  cannot  be. 

(3)  Prom  (2),  it  is  obvious  that  the  dofs  at  the  mid-side  and  center  nodes  of  elements  of  type  (a) 
and  (b)  are  different  as  well. 

(4)  A  comprehensive  treatment  of  the  mapping  procedures  for  these  elements  are  beyond  the 
scope  of  this  thesis  and  will  be  a  subject  of  future  study. 

2.4.3  Present  Study 

In  the  present  study  we  consider  C1,1  local  approximations  in  which  p- levels  in  the  £  and  y  directions 
arc  minimally  conforming  and  the  element  shapes  are  rectangular  with  the  £  and  y  axes  parallel 
to  the  x  and  y  axes.  In  this  case,  p-levels  in  both  £  and  y  arc  three.  There  are  nodal  dofs  at  the 
hierarchical  nodes.  The  dofs  at  each  of  the  corner  nodes  arc  </>,  and  We  also  consider 

p-lcvcls  beyond  minimally  conforming  in  which  case  there  are  dofs  at  the  mid-side  and  center  nodes. 
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Remarks 


(1)  Computation  of  <p ,  and  at  the  non-hierarchical  nodes  of  presents  no  problem  due  to 

global  differentiability  of  class  C1  in  x  and  y.  However,  accurate  determination 

of  requires  ml0/l  to  be  of  class  C 2  in  x  and  y.  Hence,  in  this  case  when  ml4>h  is  mapped 
onto  0^2 1  Sxdy  may  be  corrupted.  Numerical  studies  presented  in  chapter  3  demonstrate  the 
consequences  of  the  inaccuracies  in  at  the  corner  nodes  of  C^l2. 

(2)  Alternate  methodology  is  needed  to  determine  accurate  maps  of  Use  of  Lagrange  C° 

configurations  is  a  possibility. 

(3)  Numerical  studies  for  p-levels  beyond  minimally  conforming  are  also  presented  in  Chapter  3. 

2.5  Rediscretization  for  Non-linear  BVPs 

The  mapping  procedure  may  result  in  inaccurate  maps  of  the  solution  onto  rediscrctization  due 
to  either  using  non- weakly  converged  solutions  or  due  to  inaccurate  maps  of  the  hierarchical  dofs. 
This  is  critical  when  the  differential  operators  are  self-adjoint  or  non-self  adjoint  due  to  the  fact  that 
in  both  cases  the  differential  operators  are  linear  and  there  are  no  straight-forward  mechanisms  to 
improve  the  accuracy  of  the  mapped  solution.  When  the  differential  operators  are  non-linear,  the 
inaccurately  mapped  solution  can  be  used  as  a  starting  or  initial  solution  for  the  rediscrctization 
and  Newton’s  method  with  line  search  can  be  employed  to  obtain  the  improved  solution  for  the 
rediscrctization.  This  solution  so  obtained  is  expected  to  be  identical  to  the  solution  that  one 
would  obtain  if  new  computations  were  done  for  the  rediscrctization.  This  feature  is  powerful 
in  the  sense  that  it  permits  us  to  eliminate  the  mapping  errors  in  the  solution  map  onto  the 
rediscrctization.  Numerical  studies  are  presented  using  steady-state  Burgers  equation  as  a  model 
problem  to  demonstrate  this  feature  for  non-linear  BVPs. 
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2.6  Moving  Meshes  for  ID  IVPs 


In  many  IVPs,  isolated  high  solution  gradients  (fronts)  exist  in  relatively  small  portions  of  the  total 
domain  which  propagate  as  evolution  proceeds  in  time.  In  such  cases,  a  highly  refined  discretization 
may  be  necessary  in  the  portion  of  the  domain  containing  these  solution  gradients,  but  a  relatively 
coarser  mesh  may  be  employed  elsewhere.  However,  due  to  the  fact  that  the  zone  of  high  solution 
gradient  advances  in  time  during  evolution,  it  becomes  necessary  to  employ  the  same  level  of 
discretization  ahead  of  the  fronts  if  the  fronts  arc  to  be  resolved  accurately.  This  gives  rise  to 
spatial  discretization  that  may  be  computationally  impractical.  The  Riemann  shock  tube  is  a 
classic  example  of  these  types  of  problems. 

In  such  IVPs,  an  alternate  strategy  may  be  employed  to  overcome  the  problem  of  impractical 
discretizations.  The  moving  mesh  approach  is  one  such  option.  In  this  work,  we  consider  space-time 
integral  forms  of  the  IVPs  over  a  single  space-time  strip  in  which  the  space-time  integral  forms  arc 
space-time  variationally  consistent  and  employ  a  space-time  time  marching  process  to  compute  the 
evolution  state  until  the  desired  time  is  reached  [1].  The  space-time  strip  for  the  current  increment 
of  time  is  advanced  in  space  for  the  subsequent  increment  of  time  with  the  same  speed  as  that  of 
the  propagating  front  without  rediscretization.  This  process  is  continued  until  the  desired  time  is 
reached.  Details  arc  presented  in  the  following. 

Figure  2.0  shows  a  single  space-time  strip  (ID  IVP)  for  an  increment  of  time  At  for  to  <  t  < 
to  +  At,  to  being  the  initial  time  which  can  be  conveniently  assumed  to  be  zero.  Consider  the  first 
space-time  strip  discretization  (Figure  2.6  (a)).  There  are  four  distinct  zones  A-D.  Zone  B  contains 
the  initial  disturbance  defined  by  initial  conditions.  The  spatial  discretization  in  this  zone  is  such 
that  the  initial  front  is  resolved  accurately  for  a  chosen  At,  p-levcls  in  space  and  time  as  well  as 
k  —  (k\,  ka)  orders  of  the  approximation  spaces  in  space  and  time.  Zone  C  has  the  same  level  of 
refinement  as  zone  B,  but  its  length  in  space  is  determined  based  upon  the  speed  of  propagation 
of  the  front  in  zone  B  at  time  t  —  to.  A  converged  solution  of  desired  accuracy  is  computed  for 
the  first  space-time  strip.  During  the  evolution,  the  front  moves  from  zone  B  to  C  with  the  desired 
resolution  and  accuracy.  Hence,  in  zone  B  there  is  no  need  for  the  same  level  of  refinement  for 
the  next  increment  of  time.  Thus,  for  the  second  increment  of  time  the  discretization  for  the  first 
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2r„ :  open  boundary 


►k  B  -»+*-  C 


2r,:l.C. 


(b)  Second  space-time  strip 


>r4 :  open  boundary 


►k  B  ->k  C  -*k 


■r2:  IC. 


(b)  First  space-time  strip 


>r3 :  B.C. 
— *  x 


t  =  t„+  2  At 


t  =  t„+  At 


t  =  t0+2At 

'r,:  B.C. 
t  =  t0+  At 


Figure  2.6:  Moving  mesh  space-time  strips  for  first,  and  second  increments  of  time:  F  indicate  front 
location 


increment  of  time  is  employed  but  is  moved  in  space  by  a  fixed  distance  determined  by  the  speed 
of  propagation  of  the  front  in  the  x-direction.  Figure  2.6  (b)  shows  the  space-time  strip  for  the 
second  increment  of  time  (to  +  At  <  t  <  to  +  2 At). 

Computation  of  the  solution  for  the  second  space-time  strip  requires  determination  of  initial 
conditions  on  boundary  2F2  of  the  space-time  strip  from  the  solution  on  boundary  ^4  of  the  first 
space-time  strip,  i.e.  lr40/t  from  ^4  of  the  first  space-time  strip  needs  to  be  mapped  onto  2F2  of 
the  second  space-time  strip  (2r20/j).  This  problem  is  similar  to  rcdiscretization  except  that  in  this 
method  the  map  of  the  existing  solution  is  only  required  from  one  boundary  of  the  existing  space- 
time  strip  onto  another  boundary  of  the  next  space-time  strip  as  opposed  to  a  complete  map  for 
the  whole  discretization  of  the  space-time  strip.  Thus,  all  of  the  concepts  and  procedures  presented 
for  rediscretization  for  2D  BVPs  arc  applicable  here  as  well. 
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Remarks 


(1)  It  is  important  to  note  that  since  the  mapping  of  the  solution  and  its  accuracy  depends 
upon  the  nature  of  the  local  approximations  (Ct]  Lagrange,  C°  p- version  hierarchical,  higher 
order  continuity,  etc),  care  must  be  taken  to  measure  the  accuracy  of  the  solution  map  onto 
boundaries  l+1T2  from  boundaries  T4  (i  being  the  space-time  strip  number)  before  time 
marching  is  commenced.  Corrupted  solution  maps  will  obviously  result  in  faulty  evolution 
and  may  even  cause  the  computations  to  cease. 

(2)  The  thrust  of  the  work  presented  here  has  been  development  of  methodologies,  concepts  and 
their  validity.  Applications  of  these  concepts  to  complex  IVPs  will  be  the  subject  of  future 
work. 
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Chapter  3 


Numerical  Studies 


In  this  chapter  wc  present  numerical  studies  for  rcdiscrctization  for  ID  and  2D  BVPs  (self-adjoint, 
non-self  adjoint  and  non-linear  differential  operators)  and  moving  meshes  for  ID  IVP.  Numerical 
studies  are  organized  as  follows. 

(a)  Onc-dimensional  BVP  using  steady-state  ID  convection  diffusion  equation  as  a  model  prob¬ 
lem. 

(al)  Local  approximations  of  class  C°  Lagrange  type. 

(a2)  C°  p- version  hierarchical  local  approximations. 

(a3)  Local  approximations  of  type  C-7;  j  =  1,2,...  with  minimally  conforming  p-lcvels. 

(a4)  Local  approximations  of  type  C-7;  j  —  1,2,...  with  p-levels  higher  than  minimally  con¬ 
forming. 

(b)  Two-dimensional  BVP  using  Poisson’s  equation  as  a  model  problem. 

(bl)  Local  approximations  of  class  C0,0  Lagrange  type. 

(b2)  C1,1  local  approximations  with  minimally  conforming  p- level  as  well  as  p-levcls  higher 
than  minimally  conforming. 

(c)  One-dimensional  non-linear  BVP  using  ID  steady-state  Burgers  equation  as  a  model  problem. 

(d)  One-dimensional  IVP  using  ID  transient  convection  diffusion  equation. 
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3.1  One-dimensional  BVP:  Steady-state  Convection  Diffusion  Equa¬ 
tion,  Local  Approximations  of  Class  C°  Lagrange  Type 


(Pe  =  10) 

Consider 


d(f) 

dx 


1  d2<p 
Pe  dx 2 


— -  — — ^  =0  in 


0=  (0,1) 


with 

(p(0)  =  1  and  (f>(l)  —  0. 


(3.1) 


A  theoretical  solution  of  (3.1)  can  be  found  in  reference  [9],  The  differential  operator  in  (3.1)  is 
non-self  adjoint.  We  consider  a  least  squares  formulation  of  (3.1)  using  auxiliary  variable  r  =  ^  so 
that  solutions  of  class  C°  can  be  considered.  The  integral  for  in  this  case  is  variationally  consistent. 
We  consider  a  low  Pe,  say  10,  for  which  the  solution  is  highly  diffused  and  free  of  sharp  fronts  and 
hence  has  good  smoothness.  We  consider  a  ten-element  uniform  discretization  with  p-level  p  =  1 
(original  discretization:  mi).  The  rcdiscretization  considered  here  is  a  twelve-element  uniform  mesh 
(m2)  with  p- level  p  =  1.  We  consider  the  following. 

(i)  Computation  of  the  solution  for  the  original  ten-element  discretization:  TOV/t  (original  solu¬ 
tion). 

(ii)  A  map  m1^/,  onto  rcdiscretization  m2:  (mapped  solution). 

(iii)  Computation  of  a  new  solution  for  the  rcdiscretization  m2:  (m2(f)h)new  (new  solution). 

Due  to  coarser  discretization  and  low  p- level,  ml(f)h  has  significant  error.  The  least  squares  functional 
I,  dofs,  p-lcvels  and  Z/2-norms  of  the  errors  are  given  in  Table  3.1.  The  purpose  of  this  study  is 
to  show  how  the  solution  errors  map  onto  the  rcdiscretization  when  the  original  solution  is  not 
weakly  converged.  Figures  3.1  and  3.2  show  graphs  of  (f)  and  P  versus  x.  Solution  (f>  from  all  three 
cases  (original,  mapped  and  new)  arc  in  disagreement.  We  note  that  the  locations  of  the  g  jump 
discontinuities  shift  and  their  magnitudes  change  as  well.  The  mapped  solution  is  not  in  agreement 
with  the  new  solution,  which  is  the  correct  solution  for  the  rcdiscretization  m2.  The  study  clearly 
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shows  the  problems  associated  with  mapping  unconvergcd  solutions  from  the  original  discretization 
onto  the  rediscretization.  The  solution  errors  in  ml  and  m2  are  not  the  same.  The  inter-element 
flux  problems  in  ml  shift  to  different  locations  in  m2  with  different  magnitudes.  Mapping  such 


ml4>h  is  of  little  value,  if  any. 

We  now  consider  a  highly  refined  uniform  discretization  ml  of  1000  elements  with  p  —  1  and 
obtain  mV/i  and  its  map  onto  a  rediscretization  containing  1200  uniform  elements  (m2).  We  also 
compute  a  new  solution  for  the  rediscrctization  by  revolving  the  BVP.  Figures  3.3  and  3.4  show 
graphs  of  (f>  and  ^  versus  x.  In  this  case,  the  solution  is  sufficiently  converged  (weakly)  and 
thus  all  three  solutions  arc  their  derivatives  are  in  excellent  agreement.  No  visible  inter-element 
flux  problems  exist  in  any  of  the  solutions.  Table  3.2  shows  details  of  least  squares  error  functional 
I,  dofs,  error  norms,  etc. 

This  study  shows  the  importance  of  the  weak  convergence  of  ml  4>}t  before  mapping  it  onto  the 
rediscretization  m2.  Weakly  converged  solutions  when  mapped  onto  the  rediscretization  do  not 
exhibit  significant  flux  problems.  In  this  study,  p  =  1  has  been  used.  Similar  findings  also  hold  for 
higher  p-lcvcls  when  the  local  approximations  are  of  C°  Lagrange  type. 


Table  3.1:  Convection  diffusion  equation  (BVP):  C°  Lagrange  local  approximations. 


/ 

Uh  -  Hl2 

HilBllHi 

10 

1 

0.278391e0 

0.130565e0 

0.800277e0 

■m 

12 

1 

- 

0.134777e0 

0.799217e0 

{m24>h)new 

12 

1 

0.215044c0 

0.100951e0 

0.658679e0 

mL4>h 

1000 

1 

0.416684e-4 

0.196047e-4 

0.645575e-2 

1200 

1 

- 

0.202418c-4 

0.552887c-2 

(m2<Ph)neW 

1200 

1 

0.289365c-4 

0.136146c-4 

0.537967e-2 

3.2  One-dimensional  BVP:  Steady-state  Convection  Diffusion  Equa¬ 
tion,  Local  Approximations  of  C°  p-version  Hierarchical  Type 


(Pe  =  10) 

We  consider  the  BVP  (3.1).  We  consider  two  different  original  discretizations. 

Case  1:  A  tcn-elcment  uniform  original  discretization  and  a  twelve  element  uniform  rcdiscretization 
both  at  p-levels  of  5,  9,  and  15. 

Case  2:  A  1000-element  uniform  original  discretization  and  a  1200  element  uniform  rediscretization 
both  at  p- levels  of  5,  9,  and  15. 

For  both  original  discretizations,  solutions  arc  computed  at  p-levels  of  5,9  and  15  and  then  mapped 
onto  the  corresponding  rcdiscrctizations.  Also,  new  solutions  are  computed  for  the  rcdiscretizations 
by  resolving  the  BVP.  Figures  3.5  -  3.7  show  graphs  of  <fr  versus  x  for  p  =  5, 9,  and  15  for  case 
1.  The  graphs  of  ^  versus  x  are  shown  in  Fig.s  3.8  -  3.10  for  case  1.  Figures  3.5  -  3.7  show  the 
original  and  new  solutions  in  good  agreement  at  all  three  p-levels,  however  the  mapped  solution 
differs  (maybe  not  significantly).  The  graphs  of  ^  versus  x  in  Fig.s  3.8  -  3.10  clearly  show  excellent 
agreement  between  the  original  and  new  solutions  but  the  mapped  solution  differs  significantly  and 
shows  significant  inter-element  jumps  in  at  all  three  p-levels.  This  is  due  to  inaccurate  maps  of 
the  hierarchical  derivatives  which  progressively  deteriorate  as  the  p-levels  and  hence,  the  orders  of 
the  derivatives  increase.  At  all  three  p-levels  ml(j)h  is  converged,  but  the  inaccuracies  in  the  maps 
arc  due  to  inaccurate  mapping  of  the  hierarchical  dofs.  Table  3.2  gives  mesh  details,  I,  and  various 
Z/2-norms. 

The  inaccuracies  in  the  map  shown  above  can  be  minimized  with  mesh  refinement  (as  in  case 
2).  Figures  3.11  -  3.13  show  graphs  of  (f>  versus  x  and  Fig.s  3.14  -  3.16  show  graphs  of  ^  versus 
x  for  all  three  p-levels.  The  results  show  excellent  agreement  between  all  three  solutions  and  their 
derivatives  at  all  three  p-levels.  Even  though  the  L2-norms  differ  for  <p  and  ^  for  the  mapped 
solution,  for  all  practical  purposes  (as  shown  in  the  figures)  the  maps  arc  generally  good.  In  this 
case,  we  have  purposely  used  an  exceptionally  refined  mesh  to  demonstrate  that  the  inaccuracies 
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in  the  maps  of  the  hierarchical  dofs  can  indeed  be  minimized  by  mesh  refinement. 


Table  3.2:  Convection  diffusion  equation  (BVP):  C°  p- version  hierarchical  local  approximations. 


Solution 

Elements 

p-level 

/ 

II <t>h  -  4>\\l2 

d(<t>h)  d4> 

dx  dx  l 

mV/> 

10 

5 

0.432714e-9 

0.104665c-6 

0.207345c-4 

■»mr 

12 

5 

- 

0.615416e-2 

0.2360702e0 

new 

12 

5 

0.728349e-10 

0.357703e-7 

0.851514e-5 

ml4>h 

10 

9 

0.728966c-21 

0.752119c-13 

0.269668el0 

mim 

12 

9 

- 

0.615716e-2 

0.236839e0 

fttfrrwni 

12 

9 

0.286326c-22 

0.125148e-13 

0.534643e-ll 

"i '  ™ii 

15 

0.944106c-22 

0.252814e-13 

0.914020e-ll 

mmm 

12 

15 

0.615718e-2 

0.409828e-2 

bbsim 

12 

15 

0.205984c-21 

0.299340c-13 

0.132347c-14 

HIT 

1000 

5 

0.296098c-23 

0.400021e-12 

0.367759e-ll 

fT 

1200 

5 

- 

0.623740e-6 

0.236693e-2 

bbimi 

1200 

5 

0.551350e-21 

0.109110c-10 

0.443337e-10 

wumm 

1000 

9 

0.6334062e21 

0.106185c-10 

0.463778e-10 

— 

9 

- 

0.623748c-6 

0.236693e-2 

(m^<fih)new 

1200 

9 

0.926781e-21 

0.139454e-10 

0.584280e-10 

m 

1000 

15 

0.247810e-17 

0.147141e-10 

0.156245c-8 

1200 

15 

- 

0.789122e-6 

0.409828e-2 

(m‘2<Pk)new 

1200 

15 

0.141826e-17 

0.297155e-ll 

0.116684e-8 
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Original 

mesh 

0 

C  ,  p=9. 

10  el . 

Rediscretization 

c°,  p=9, 

12  el . 

ml  , 

4*h  ' 

m2  , 

<i>h ' 

.m2 1  . 

'  ■  h '  new 

Figure  3.6:  Convection  diffusion  equation  (BVP):  0  versus  x  for  original,  mapped  and  new  solutions. 


Solution  d(|>/dx 


Figure  3.8:  Convection  diffusion  equation  (BVP):  ^  versus  x  for  original,  mapped  and  new  solu¬ 
tions. 
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d(j)/dx 


d0/dx 


3.3  One-dimensional  BVP:  Steady-state  Convection  Diffusion  Equa¬ 
tion,  Local  Approximations  of  0\  j  =  1,2, 3, 4  ( Pe  =  10) 

In  this  case  we  consider  the  original  discretization  to  be  a  five-element  uniform  mesh  (ml).  The  rc- 
discretization  consists  of  s  seven-element  uniform  mesh  (m2).  In  the  studies,  we  consider  minimally 
conforming  p-levels  for  each  class  as  well  as  p- levels  higher  than  minimally  conforming.  Results  are 
presented  in  the  following  figures  and  tables. 

(a)  C1,  p  =  3:  Figures  3.17  -  3.19,  Table  3.3 

(b)  C1,  p  =  15:  Figures  3.20  -  3.22,  Table  3.3 

(c)  C2,  p  —  5:  Figures  3.23  -  3.26,  Table  3.4 

(d)  C2,  p  =  15:  Figures  3.27  -  3.30,  Table  3.4 

(c)  C3,  p  =  7:  Figure  3.31,  Table  3.5 

(f)  C3,  p  =  15:  Figure  3.32,  Table  3.5 

(g)  C4,  p  =  9:  Figure  3.33,  Table  3.6 

(h)  C4,  p  =  15:  Figure  3.34,  Table  3.6 

First,  we  discuss  results  for  with  minimally  conforming  p-levcls.  The  solution  of  progressively 
higher  classes  exhibit  progressively  better  maps  of  the  solution  as  well  as  the  derivatives  of  up  to 
order  C3 .  When  the  solutions  arc  of  class  C3 ,  the  maps  of  the  derivatives  of  order  j  +  1  improve 
with  progressively  increasing  order  of  the  approximation  space  when  the  p-levels  are  minimally 
conforming.  This  is  due  to  increasing  p-lcvcls  progressively  by  increasing  the  order  of  the  space  but 
with  the  absence  of  hierarchical  dofs  which  may  not  map  accurately.  When  p-levcls  are  increased 
beyond  minimally  conforming,  the  map  of  the  derivatives  of  order  C3+1  become  inaccurate  (Fig.s 
3.22,  3.30)  for  classes  C 1  and  C2  but  remain  good  for  classes  C3  and  C4.  This  is  due  to  the  fact 
that  minimally  conforming  p-levels  for  classes  C3  and  C4  are  7  and  9,  respectively,  at  which  the 
original  solution  is  relatively  well  converged  and  there  arc  still  no  hierarchical  dofs.  Hence,  the 
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maps  onto  the  rediscretization  should  be  good.  At  p-level  15,  ^  and  show  very  good  maps 
for  solutions  of  classes  C3  and  C"4  (Fig.s  3.32,  3.34). 

Table  3.3:  Convection  diffusion  equation  (BVP):  C 1  p-version  hierarchical  local  approximations. 


Solution 

Elements 

p-level 

I 

\\<t>h  -  4>\\L2 

d(4>h)  d<ji 

dx  dx  j 

<P(4>h)  d?<t> 

dx 2  dx^  ^ 

5 

3 

0.730018c-l 

0.295100e-l 

0.132388e0 

0.321291el 

~^h 

7 

3 

- 

0.290958e-l 

0.126901e0 

0.307579el 

(m‘2j>h)neW 

7 

3 

0. 234047c- 1 

0.951054c-2 

0.465439e-l 

0.131715el 

5 

15 

0.847645e-23 

0.175931e-14 

0.180028e-l 

0.290602e-10 

■  gj  ■ 

7 

15 

- 

0.406817e-2 

0.988658e-l 

0.538991el 

MSBENMM 

7 

15 

0.382862c-23 

0.323930e-14 

0.892078e-13 

0.195499e-10 
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Table  3.5:  Convection  diffusion  equation  (BVP):  C3  /^version  hierarchical  local  approximations. 


Solution 

Elements 

p- level 

I 

d4(<Ph)  dU 

dx 4  dx4  £ 

WBBEBM 

5 

7 

0.448956c-8 

0.976222el 

bbhw 

7 

7 

- 

0.937789el 

WBmrnm 

7 

7 

0.115551e-9 

0.281397el 

Banff 

5 

15 

0.129770c-22 

0.606670e-5 

«a mm 

7 

- 

0.228179e2 

C m'2<t>h)nem 

7 

15 

0.821314e-22 

0.261075e-4 

Table  3.6:  Convection  diffusion  equation  (BVP):  CA  p- version  hierarchical  local  approximations. 


Solution 

Elements 

p-levcl 

I 

■ 

ml4>h 

5 

9 

0.153421e2 

— 

7 

9 

- 

0.148577e2 

(’ m24>h)new 

7 

9 

0.608197e-18 

0.310044el 

ml4>h 

5 

15 

0.107051e-22 

0.176433e-2 

m2<»h 

7 

15 

- 

0.228179e2 

(m‘2MneW 

7 

15 

0.260707e-23 

0.161274e-2 
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Ill 


d  0/dx 


xp/<t)cp 


x 

Figure  3.34:  Convection  diffusion  equation  (BVP):  versus  x  for  original,  mapped  and  1 

solutions. 
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3.4  Two-dimensional  BVP:  Poisson’s  Equation,  Local  Approxima¬ 
tions  of  Classes  C0  0  Lagrange  Type  with  p^  =  pv  =  1,15  and 
C1,1  with  =  p^  =  3, 15 

Wc  consider  steady-state  Poisson’s  equation 

<™er  fi  =  (-l.1)*(-l»l)  (3-2) 

with 

0(-i,  y)  =  0(1,  y)  =  <t>(x>  -1)  =  0(*,  i) 

where  f(x,y )  is  such  that  the  theoretical  solution  is  given  by 

(p{x,y)  =  e^1_a;2hi-y2). 

The  differential  operator  is  self-adjoint  and  hence  the  integral  form  from  Galerkin  method  with 
weak  form  is  variationally  consistent.  Figure  3.35  (a)  shows  a  schematic  of  the  domain  fl  as  the 
original  5x5  uniform  discretization  (Figure  3.35  (b))  as  well  as  a  7  x  7  uniform  rcdiscretization 
(Figure  3.35  (c)).  Wc  consider  the  following  studies. 

(1)  Local  approximations  of  class  C0,0  Lagrange  type  with  =  pv  =  1  and  of  class  C°'°  p- version 
hierarchical  with  p-levcl  of  p^  =  p^  =  15.  Wc  consider  the  integral  form  based  on  Galerkin 
method  with  weak  form. 

(2)  Local  approximations  of  class  C1,1  p- version  hierarchical  with  p-levcls  of  =  p,}  =  3, 15.  The 
integral  form  is  based  on  least  squares  method  without  auxiliary  variables  (i.e.  strong  form 
of  the  GDE). 

In  all  cases  ml0fc,  is  computed  for  ml  and  mapped  onto  m2  to  obtain  m20/(.  A  new  solution  is  also 
computed  for  m2  to  obtain  ( m2<ph)new  for  the  rediscretization  m2  by  resolving  the  BVP.  Results 
are  summarized  in  the  following  figures  and  tables. 

(a)  G0,0,  p(_  —  pv  —  1;  Figures  3.36,  3.37  (a),  (b)  and  (c),  Table  3.7 
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(a)  Schematic  of  Q  for  Poisson’s  equation 


(b)  5x5  uniform  discretization  ml  (c)  ^x7  uniform  rediscretization  m2 

Figure  3.35:  Schematic,  discretization  and  rcdiscretization  for  Poisson’s  equation. 
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(b)  C0,0,  p£  —  pv  —  15;  Figures  3.37  (d),  (e)  and  (f),  Table  3.7 

(c)  C1,1,  p^  —  pv  —  3;  Figures  3.38,  3.39  and  3.40,  Table  3.8 

(d)  C'1,1,  p^  —  Prj  =  15;  Figures  3.41,  3.42  and  3.43,  Table  3.8 

The  ml(f>h  solution  of  class  C0,Q  with  =  pv  —  1  for  discretization  ml  is  unconverged  and  hence 
has  errors  (Table  3.7).  Graphs  of  ^  and  clearly  demonstrate  inter-element  jumps  in  the  first 
partial  derivatives  at  different  locations  for  and  m20/,  with  different  magnitudes  as  well. 
is  not  in  agreement  with  ( ml(ph)newi  which  is  the  correct  solution  for  m2.  This  demonstrates  the 
inter-element  flux  problems  caused  in  the  mapped  solution  due  to  inter-element  flux  problems  in 
the  original  solution  ml^l.  Due  to  p ^  —  p7)  =  1  only  function  values  are  dofs  at  the  nodes.  Just 
as  in  the  case  of  the  ID  convection  diffusion  equation,  here  also,  mesh  refinement  will  minimize 
the  flux  problems  in  ml4>h  as  well  as  m2<ph-  Higher  p- level  Lagrange  elements  behave  in  a  similar 
fashion  similar  to  the  ID  case. 

Next,  wc  consider  C°  p- version  hierarchical  approximations  for  ml  and  m2  (Figures  3.37  (d)- 
(f)).  Even  though  the  mesh  is  relatively  coarse,  due  to  high  p-level  ml4>h  and  ( m24>h)new  arc 
sufficiently  converged  and  are  in  good  agreement  with  each  other  (Figures  3.37  (d)  and  (f)).  Due 
to  the  C°  nature  of  the  local  approximations,  maps  of  the  hierarchical  dofs  at  the  mid-side  and 
center  nodes  of  the  elements  require  maps  of  second  and  higher  order  derivatives  of  <p  with  respect 
to  £  and  p.  As  a  result,  m2<Ph  is  spurious  (Figure  3.37  (c)).  Similar  behavior  is  also  observed 
for  the  maps  of  the  derivatives  and  ^  onto  rediscretization  m2  (not  shown).  This  problem 
can  be  minimized  or  possibly  corrected  by  first  mapping  the  solution  onto  the  corresponding 
Lagrange  configurations  (requiring  only  function  values)  and  then  obtaining  the  hierarchical  dofs 
from  these  configuration. 

The  solutions  of  class  C1,1  at  =  p^  =  3  (Figures  3.38-3.40)  show  good  maps  and  are  in  good 
agreement  with  ml0/i  and  ( ml<frh)new ,  indicating  no  potential  mapping  problems  for  cj)  as  well  as  its 
derivatives  up  to  second  order  confirming  weak  convergence  of  the  original  solution  rnl(ph  as  well 
as  ( ml(ph)new  Solutions  of  class  C1,1  at  =  pv  =  15  (Figures  3.41-3.43)  show  good  maps  as  well 
even  though  higher  order  hierarchical  dofs  arc  mapped  in  this  case.  This  is  due  to  the  fact  that 
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since  =  pn  =  3  solutions  (with  no  hierarchical  dots)  are  sufficiently  converged,  the  hierarchical 
dofs  beyond  p ^  =  pn  =  3  are  not  contributing  significantly  to  the  improvement  in  the  solution  and 
hence  the  inaccuracies  in  their  maps  arc  of  little  or  no  consequence. 


Table  3.7:  Poisson’s  equation  (BVP):  C0,0  Lagrange  local  approximations. 


Solution 

Mesh 

p- level 

I 

Hollis 

IBM 

1— 

5x5 

1 

- 

0.336199cl 

0.247853el 

wmm- 

1 

- 

0.331366cl 

0.237681el 

(m'20h)nei« 

7x7 

1 

- 

0.338225el 

0.250033cl 

wkmmh 

- 

0.3403821el 

0.252413el 

fmrm 

15 

- 

0.375170el 

0.105981e2 

15 

- 

0. 340557c! 

0.252327el 

Table  3.8:  Poisson’s  equation  (BVP):  C1,1  p-version  hierarchical  local  approximations. 


Solution 

I 

\\<t>h\\L2 

,am 

d2{<t>h) 
dx'2  L, 

dxdy  Lf> 

iny 

5x5 

3 

0.368110c-l 

0.340376cl 

0.252341el 

0.424995cl 

0.401647el 

TW* 

7x7 

3 

- 

0.340365cl 

0.252322el 

0.424924cl 

0.401610el 

(m2<t>k)neW 

7x7 

3 

0.964658c-2 

0.340380el 

0.252394cl 

0.425025el 

0.401784cl 

mBEm 

5x5 

15 

0.321561c-ll 

0.340382cl 

0.252413el 

0.425034el 

0.401835cl 

m 

15 

- 

0.340363cl 

0.252359el 

0.425034el 

0.401708el 

(m^4>h)new 

7x7 

15 

0.100513e-9 

0. 340384c  1 

0.252411el 

0.425035el 

0.401847el 
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3.5  One-dimensional  Non-linear  BVP:  Steady-state  Burgers  Equa¬ 


tion  (Re  —  25) 


Here,  wc  consider  the  steady-state  Burgers  equation  as  a  model  BVP. 


d<!> 


1  d2(j) 
Pe  dx 2 


=  0  in 


O  =  (0,1) 


(3.3) 


with 

0(0)  =  1  ,  0(1)  =  0. 

We  consider  Re  =  25,  for  which  the  theoretical  solution  is  smooth  and  there  arc  no  sharp  fronts. 
We  consider  a  uniform  discretization  of  five  p-version  elements  (ml)  and  a  uniform  rediscrctization 
of  seven  p- version  elements  (m2).  In  (3.3),  the  differential  operator  is  non-linear.  We  consider 
least  squares  formulation  of  (3.3)  (strong  form  of  GDE)  for  which  C 1  is  the  minimally  conforming 
approximation  class  if  wc  permit  weak  convergence  of  the  second  derivative  of  0.  We  consider  the 
following  two  cases. 

(a)  p-levcl  of  5,  at  which  ml0/t  is  unconvcrged  and  hence  has  errors. 

(b)  p-levcl  of  15,  at  which  the  solution  ml0/,  is  weakly  converged  and  hence  is  relatively  free  of 
errors. 

For  each  of  the  above  cases,  wc  choose  local  approximations  of  class  C 1  and  perform  the  following 
computations. 

(i)  A  solution  is  computed  for  discretization  ml:  ml0/t- 

(ii)  ml0/»  is  mapped  onto  to  obtain  m20/t. 

(iii)  A  new  solution  is  computed  for  the  rediscretization  m2  using  boundary  conditions  in  (3.2): 

(m20/>W 

(iv)  The  mapped  solution  ra20/l  is  used  as  an  initial  solution  for  rediscretization  m2  and  Newton’s 
method  with  line  search  is  employed  to  obtain  a  converged  solution:  {m24>h)new  by  resolving 
the  BVP. 
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The  results  are  summarized  in  the  following  figures  and  Table  3.9. 

(a)  C\  ^  =  3;ml4,m24,(ml4)  new  and  their  derivatives  up  to  second  order;  Figures  3.44-3.46 

(b)  C\  Pi  =  3  ;  m24>lu  ("*Vfc)  new j  (ml(t>h)new  and  their  derivatives  up  to  second  order;  Figures 
3.47-3.49 

(c)  C\  =  15  ; 

new  and  their  derivatives  up  to  second  order;  Figures 

3.50-3.52 

(d)  C1,  p£  —  15  ;  (ml4>h)new,  (,ml<l>h)new  and  their  derivatives  up  to  second  order;  Figures 

3.53-3.55 

Since  at  —  3  the  solution  m  1  is  not  converged,  it  has  errors  and  hence  its  map  rn2(ph  is 
expected  to  have  errors  as  well  and  the  maps  of  the  errors  in  the  two  cases  will  also  be  different. 
In  Figure  3.44,  the  solution  map  is  in  fair  agreement  with  ml<f>h ,  however  the  correct  solution  for 
m2  is  (m20/i)netti  which  does  not  agree  with  the  map.  The  differences  between  the  mapped  solution 
rn2<ph  and  (m2</>/l)ne ui  are  more  pronounced  in  Figures  3.45  and  3.46  showing  graphs  of  the  first 
and  second  derivatives,  respectively.  From  Figures  3.47-3.49,  showing  graphs  of  m2^,  (m2q bfl)new, 
and  (ml (fih) new  and  their  derivatives  up  to  second  order,  we  observe  perfect  agreement  between 
(  < Ph)new  and  and  their  derivatives  up  to  second  order.  Even  though  rn2(ph  is  in  error, 

using  m2<^)/i  as  an  initial  solution,  it  is  possible  to  obtain  the  correct  solution  for  the  rediscretization 
m2.  This  is  a  significant  merit  of  the  proposed  approach  using  as  an  initial  solution  for  the 
BVPs  described  by  non-linear  differential  operators.  Solutions  of  class  Cl  at  pi  =  15  (Figures 
3.50-3.55)  further  substantiate  the  usefulness  and  merit  of  this  approach.  From  these  figures,  we 
note  that  even  though  m2<f>h  and  their  derivatives  are  erroneous,  (m2 tyhYnew  in  perfect  agreement 
with  (m20/i)ne w  and  their  derivatives  of  up  to  second  order.  The  findings  reported  here  also  hold 
when  the  local  approximations  are  of  higher  classes. 
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Tabic  3.9:  Burgers  equation  (BVP):  C1,1  p-version  hierarchical  local  approximations. 


p- level 

I 

3“ 

to 

1 

IHsUll  SI 

1 

ml<t>h 

5 

3 

0.105577e0 

0.829991e0 

o 

.183169el 

0.131250e2 

wmm 

7 

3 

- 

0.830241e0 

0.183269el 

0.130113c2 

7 

3 

0.799813c-l 

0.850746c0 

0.199550el 

0.155220e2 

7 

3 

0. 79981 3c- 1 

0.850746c0 

0.199549el 

0.155220e2 

5 

15 

0.113703e-10 

0.959166e0 

0.288675el 

0.322748e2 

MSM2SM 

7 

15 

- 

0.963128c0 

0.317969el 

0.495645c2 

(m'Z<t>h)  new 

7 

15 

0.154861e-12 

0.959166c0 

0.288675el 

0.322748e2 

Kf 

7 

15 

0.155708c-12 

0.959166e0 

0. 288675c! 

0.322748c2 
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Figure  3.44:  Burgers  equation  (BVP):  <j>  versus  x  for  original,  mapped  and  new  solutions. 


Original  mesh 
C1,  p=15,  5  el. 

Rediscretization 
C1,  p=15,  7  el. 


(m2<i>h)Tiew,  r^h) 


m2 i  v is 


h'  nev/ 


h'  new 


3.6  Moving  Mesh  ID  IVP:  Transient  Convection  Diffusion  Equa¬ 


tion  ( Pe  =  106) 


We  consider  the  following  transient  convection  diffusion  equation 


d(p  d(j)  1  d2(j)  . 

fa  +  -di~p^W=  m 


^  x  0*  =  (0, 1)  x  (0,  t) 


(3.4) 


with 


0(O,f)  =  O<?!>'(O,f)  =  O 

4>(x,0)  =  exp(-~-^~)  (3.5) 

/<T0 

where  Xq  and  ctq  arc  the  mean  and  standard  deviation  of  the  Gaussian  distribution  and  arc  chosen 
to  be  0.2  and  0.03,  respectively.  We  consider  the  first  and  subsequent  space-time  strips  shown  in 
Figure  3.56.  Discretization  details  for  the  first  and  subsequent  space-time  strips  are  also  shown  in 
Figure  3.56.  The  spacial  discretization  in  zone  B  is  chosen  such  that  the  initial  conditions  described 
by  the  Gaussian  distribution  arc  resolved  accurately.  Upon  evolution,  this  disturbance  (I.C.)  moves 
into  zone  C  in  which  the  discretization  is  of  the  same  refinement  as  in  zone  B  and  hence  will  yield 
time  accurate  evolution.  A  time  step  At  of  0.02  is  considered.  For  the  second  space-time  strip  the 
discretization  is  moved  in  the  spatial  direction  by  Ax  —  0.02,  which  is  in  conformity  with  the  speed 
of  propagation  of  the  disturbance.  The  solution  from  the  open  boundary  of  the  first  space-time 
strip  is  mapped  as  I.C.  for  the  second  space-time  strip  and  the  time  evolution  is  computed.  This 
process  is  continued  for  each  increment  of  time  until  t  =  t  is  reached  (r  =  6Af  =  0.12).  Since 
the  local  approximations  arc  of  class  C1,1  with  p-level  of  three  in  space  and  time,  there  arc  no 
hierarchical  dofs  at  the  mid-side  and  center  nodes  of  the  space-time  elements,  and  furthermore, 
since  the  solution  for  the  first  space-time  strip  is  converged,  we  expect  the  I.C.  map  of  the  solution 
for  all  space-time  strips  to  be  good.  Time  evolution  of  the  solutions  is  shown  in  Figure  3.57.  At 
Pe  =  106,  the  physical  diffusion  is  insignificant  and  hence  we  expect  the  Gaussian  distribution  to 
march  in  time  without  amplitude  decay  or  base  elongation  as  shown  in  Figure  3.57. 

Remarks 
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(a)  The  moving  mesh  procedure  permits  relatively  coarser  discretization  for  Slx  of  flxy  =  flxxQt- 

(b)  Highly  refined  meshes  in  the  area  of  high  solution  gradients  and  relatively  coarser  meshes 
elsewhere  are  more  practical  and  are  in  conformity  with  what  is  necessitated  by  the  physics. 

(c)  In  this  approach,  a  micro  disturbance  and  its  propagation  can  be  resolved  over  a  macro 
domain  without  excessive  computational  effort. 

(d)  The  approach  requires  accurate  mapping  of  the  I.C.  regardless  of  the  nature  of  the  differential 
operator.  The  initial  solution  approach  presented  for  non-linear  BVPs  to  improve  the  mapped 
solution  cannot  be  used  for  I  VPs  and  moving  mesh  approach  due  to  the  fact  that  I.C.  arc 
fixed  and  hence  cannot  be  altered  during  evolution.  This  does  present  a  problem  when  p- 
levels  higher  than  minimally  conforming  for  the  chosen  order  of  the  approximation  space  are 
used  due  to  inaccuracies  in  the  map  of  the  hierarchical  dofs.  Further  work  to  alleviate  this 
problem  is  in  progress.  Initial  studies  are  promising. 
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3  unif.  el.  25  unif.  el. 


3  unif.  el. 


Figure  3.56:  First  three  moving  mesh  space-time  strips  for  transient  convection  diffusion  equation. 


9' 
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igure  3.57:  Iransient  convection  dinusion  equation 


Chapter  4 


Summary  and  Conclusions 


Rcdiscretization,  moving  mesh  and  solution  mapping  approaches  arc  presented  for  BVPs  and  IVPs. 
These  methodologies  essentially  require  two  basic  steps.  In  the  first  step,  the  geometry  of  the  new 
discretization  is  mapped  onto  the  existing  discretization  .  In  the  second  step,  the  solution  from  the 
existing  discretization  is  mapped  onto  the  new  discretization  or  rediscrctization  using  the  geometry 
map  form  step  one.  In  the  following,  we  summarize  the  work  presented  here  and  draw  some 
conclusions. 

(1)  In  the  case  of  BVPs,  rcdiscretization  requires  a  total  map  of  the  new  discretization  onto  the 
existing  mesh  and  then  a  total  map  of  the  solution  from  the  current  or  existing  discretization 
to  the  new  discretization  or  rcdiscretization. 

(2)  In  the  case  of  IVPs,  using  space-time  strips  and  space-time  time  marching,  only  a  map  of 
the  geometry  of  the  boundary  requiring  initial  conditions  onto  the  the  open  boundary  of  the 
previous  or  existing  space-time  strip  is  required.  Using  this  map,  the  computed  solution  at 
the  open  boundary  of  the  current  space-time  strip  is  mapped  onto  the  boundary  of  the  moved 
space-time  strip  that  requires  initial  conditions. 

(3)  Except  the  basic  difference  discussed  in  (1)  and  (2),  the  mapping  procedures  for  the  geometry 
and  the  solution  for  rcdiscretization  (for  BVPs)  and  moving  meshes  (IVPs)  are  the  same. 

(4)  First  we  consider  rcdiscretization  for  BVPs. 
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(a)  Mapping  of  the  rediscretization  geometry  onto  the  existing  discretization  is  always  exact 
and  unique  when  quadratic  mapping  is  assumed  for  the  subdomains  between  the  physical 
coordinate  space  x,  ( x ,  y)  or  ( x ,  y,  z )  and  £,  (£,  rj)  or  (£,  77,  Q 

(b)  Accuracy  of  the  map  of  the  solution  from  the  existing  discretization  onto  the  rediscretiza¬ 
tion  depends  upon  the  type  of  local  approximations. 

(bl)  When  local  approximations  are  of  class  C°  Lagrange  type  in  which  only  function 
values  arc  nodal  degrees  of  freedom,  the  solution  mapping  procedure  presented  here 
is  good  and  free  of  any  assumptions  or  approximations.  If  an  unconverged  solution 
is  mapped,  then  the  error  maps  in  the  two  discretizations  arc  not  unique.  More 
specifically,  in  the  case  of  C°  local  approximations,  magnitudes  of  inter-element 
jumps  in  the  solution  derivatives  normal  to  the  inter-element  boundaries  and  their 
locations  in  the  two  maps  are  different.  Thus,  inter-element  flux  problems  exist 
in  both  solutions  as  widely  reported  in  the  literature.  We  clearly  observe  that 
inter-element  flux  problems  in  the  mapped  solution  arc  to  to  the  inter-element  flux 
problems  in  the  solution  for  the  original  discretization.  When  weakly  converged 
solutions  are  mapped,  maps  of  the  solution  and  its  first  derivative  remain  good, 
implying  a  lack  of  serious  inter-element  flux  problems.  These  findings  hold  true  for 
p  =  1  as  well  as  higher  p-level  Lagrange  elements  in  ID  as  well  as  2D. 

(b2)  When  the  local  approximations  are  of  class  C°  p-version  hierarchical,  the  solution 
map  onto  the  rediscrctizations  require  mapping  of  higher  order  derivatives  with 
respect  to  £  and  p  (ID  and  2D  BVPs).  Due  to  the  C°  nature  of  the  global  ap¬ 
proximations,  direct  mapping  of  the  degrees  of  freedom  using  local  approximations 
of  the  subdomains  of  the  original  discretization  leads  to  substantial  inaccuracies  in 
the  mapped  solution.  These  inaccuracies  or  errors  increase  with  increasing  p-lcvcl 
(corresponding  to  higher  order  derivatives  with  respect  to  £  and  77).  An  alternative 
procedure  has  been  proposed  to  overcome  this  difficulty,  in  which  the  solution  is 
first  mapped  onto  the  corresponding  C°  Lagrange  configurations  that  only  requires 
function  values  and  then  the  hierarchical  dofs  are  recovered  from  these  configura- 
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tions.  With  this  procedure,  we  expect  solution  maps  for  C°  p- version  hierarchical 
local  approximations  to  be  good. 

(b3)  In  the  case  of  i  —  1,2,...  higher  order  continuity  ID  local  approximations, 
there  arc  no  problems  with  the  solution  map  when  minimally  conforming  p-levels 
are  used,  p-levels  higher  than  minimally  conforming  yield  inaccurate  solution  maps 
due  to  inaccuracies  in  the  mapping  of  hierarchical  dofs.  This  problem  can  also  be 
corrected  using  the  corresponding  C°  Lagrange  configurations  (work  in  progress). 
In  the  case  of  2D  higher-order  continuity  local  approximations,  preliminary  studies 
indicate  that  for  C1'1  solutions,  the  mapping  inconsistencies  in  ^  at  the  corner 
nodes  are  minimized  when  minimally  conforming  p-levels  are  used.  However,  for  C2,2 
and  higher  order  continuity  approximations,  these  inconsistencies  occur  in  more  dofs 
at  the  corner  nodes  even  when  minimally  conforming  p-levels  are  used.  Further  work 
is  needed  to  resolve  these  issues  (work  in  progress). 

(5)  In  the  case  of  moving  meshes,  the  computations  work  exceptionally  well  when 

(i)  Only  converged  solutions  from  the  current  space  time  strip  are  used  to  determine  ini¬ 
tial  conditions  for  the  subsequent  space-time  strip  in  which  the  mesh  has  been  moved 
spatially. 

(ii)  The  solution  mapping  process  is  free  of  errors  and  inaccuracies. 

(iii)  The  initial  conditions  are  resolved  accurately  and  sufficient  resolution  is  provided  ahead 
of  the  front  for  its  accurate  time  evolution  over  the  space-time  strip. 

(6)  When  the  mapping  procedure  for  BVPs  results  in  errors  in  the  mapped  solution  caused  by 
using  non-wcakly  converged  solutions  or  due  to  inaccurate  maps  of  the  hierarchical  dofs,  the 
situation  is  critical  when  the  differential  operators  are  self-adjoint  or  non-self  adjoint  due  to 
the  fact  that  in  both  bases  the  differential  operators  are  linear  and  there  is  no  further  mech¬ 
anism  to  improve  the  mapped  solution.  When  the  differential  operators  are  non-linear,  the 
inaccurate  mapped  solution  can  be  used  as  a  starting  or  initial  solution  for  the  rediscretiza¬ 
tion  and  the  Newton’s  method  with  line  search  can  be  used  to  obtain  the  new  or  ’correct’ 
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solution  for  the  rediscretization.  This  solution  always  identical  to  the  solution  that  one  would 
obtain  if  new  computations  were  done  for  the  rediscretization.  This  feature  is  powerful  in  the 
sense  that  it  permits  us  to  eliminate  the  mapping  errors  in  the  solution  map  onto  the  redis¬ 
cretization.  In  the  case  of  IVPs,  the  situation  is  slightly  different  due  to  the  fact  that  there 
appears  to  be  no  such  mechanism  to  eliminate  the  inaccuracies  in  the  map  of  the  hierarchical 
dofs  for  the  initial  conditions  due  to  the  fact  that  initial  conditions  are  fixed.  Thus,  for  IVPs 
an  accurate  map  of  the  initial  conditions  for  the  moved  mesh  is  essential  for  accurate  time 
evolutions. 

(7)  Benefits  of  higher-order  continuity  local  approximations  are  clearly  demonstrated.  When 
local  approximations  of  class  CJ  arc  weakly  converged,  the  solution  maps  contain  accurate 
values  of  the  solution  and  its  derivatives  up  to  order  j  +  1 .  This  is  only  possible  in  the  h-p-k 
framework. 

(8)  Numerical  studies  presented  for  ID  and  2D  BVPs  and  ID  IVPs  demonstrate  various  features  of 
rediscrctizations  and  moving  mesh  procedures  proposed  in  this  work.  Possible  remedies  have 
been  suggested  to  overcome  the  mapping  inaccuracies  in  hierarchical  dofs  and  higher-order 
continuity  local  approximations.  Preliminary  studies  arc  promising.  With  the  procedures 
proposed  here,  large  deformation  large  strain  studies  for  solid  mechanics  can  be  done  accu¬ 
rately  using  the  total  Lagrangian  approach.  Moving  fronts  for  IVPs  can  be  simulated  using 
coarser  spatial  discretizations.  Lastly,  it  is  envisioned  that  simulations  of  moving  micro  fronts 
on  a  macro  scale  will  be  possible  using  the  moving  mesh  procedure  proposed  here.  Simulation 
of  moving  shocks  in  a  Riemann  shock  tube  of  actual  physical  dimensions  is  an  example  of 
such  physics. 
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